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Use  of  measurement  models  is  suggested  as  an  intermediate  step  in  the  analysis  of 
impedance  spectra.  A measurement  model,  based  on  Voigt  circuit  elements,  is  shown  to 
apply  to  a wide  variety  of  typical  electrochemical  impedance  spectra.  The  characterization 
of  measurement  errors,  through  the  use  of  a measurement  model,  is  shown  to  help  guide 
design  of  experiments  and  model  identification. 

A new  algorithm  is  outlined  which  uses  the  measurement  concept  to  check  for 
consistency  of  time  dependent  data  with  the  Kramers-Kronig  relations.  A technique  to 
determine  stochastic  errors  for  nonstationary  systems  is  presented.  An  empirical  model  for 
the  error  structure  of  impedance  data  is  proposed.  Software,  developed  in-house,  is 
described  which  can  be  used  to  collect  impedance  data.  The  unique  features  of  the 
software  are  the  automatic  selection  of  measuring  resistor  and  an  algorithm  to  improve 
impedance  experiments  in  galvanostatic  mode. 
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The  applicability  of  electrochemical  impedance  spectroscopy  as  a viable  technique 
to  study  metal  hydride  electrodes  is  demonstrated.  The  measurement  model  approach  can 
also  be  applied  to  other  generalized  transfer  functions.  The  measurement  models  are  used 
characterize  electro-hydrodynamical  impedance  measurements. 
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CHAPTER  1 


INTRODUCTION 

Impedance  Spectroscopy  has  been  shown  to  be  a powerful  technique  for  the 
characterization  of  electrochemical  systems  [1-3],  The  promise  of  impedance  spectroscopy 
is  that,  with  a single  experimental  procedure  encompassing  a sufficiently  broad  range  of 
frequency,  the  influence  of  the  governing  physical  and  chemical  phenomena  may  be  isolated 
and  distinguished  at  a given  applied  potential  [4],  However,  like  any  other  electrochemical 
technique,  impedance  spectroscopy  cannot  be  regarded  as  being  a “stand-alone”  method 
for  model  identification  because  the  response  of  a system  to  a periodic  perturbation  does 
not  provide  a direct  measure  of  the  governing  physical  phenomena.  Interpretation  of 
impedance  data  requires  identification  of  appropriate  models.  This  inferential  character  of 
impedance  spectroscopy  and  the  ambiguity  inherent  in  the  interpretation  of  impedance 
spectra  give  rise  to  several  interesting  questions  during  model  identification.  For  example: 

•Do  the  data  satisfy  the  Kramers-Kronig  transforms  [5,6]  (i.e.,  are  the  data 
representative  of  a stationary,  causal,  and  linear  system)? 

•What  is  the  “correct”  model  for  a data  set  (i.e.,  should  a data  set  be  interpreted  in 
terms  of  coupled  reactions,  surface  roughening,  frequency  dispersion,  or  something  else)? 

•Do  the  models  really  provide  a statistically  significant  fit  to  the  data? 
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•What  portion  of  the  spectrum  can  be  interpreted  in  terms  of  a model  that  assumes 
stationary  behavior? 

•Is  the  experimental  protocol  optimized  with  respect  to  minimizing  errors  in  the  data? 

The  object  of  this  thesis  is  to  provide  a new  framework  for  model  identification  which 
emphasizes:  use  of  deterministic  models  based  on  solution  of  the  equations  governing  the 
physics  and  chemistry  of  a given  system,  application  of  experimental  techniques  sufficient  to 
confirm  hypothesized  mechanisms,  and  application  of  statistical  tools  to  ensure  that 
measurement  characteristics  can  be  taken  into  account  when  comparing  models  to 
experimental  data. 

1.1  A Framework  for  Model  Identification 

The  traditional  approach  of  characterization  of  electrochemical  system  is  illustrated  in 
the  flow  diagram  shown  in  Figure  1.1  [7],  The  experimental  data  obtained  for  the  material- 
electrode  system  is  analyzed  by  using  an  exact  mathematical  model  based  on  plausible  physical 
theory  that  predicts  theoretical  impedance  or  by  an  equivalent  circuit.  The  parameters  of  the 
model  are  obtained  by  fitting  the  model  to  the  experimental  data  by  some  regression  algorithm, 
for  example  CNLS  [8], 

The  framework  for  this  thesis  is  illustrated  schematically  in  Figure  1.2.  In  this  work  the 
general  approach  to  impedance  spectroscopy  is  divided  into  three  main  categories: 

1.  Optimal  design  of  experiments 

2.  Characterization  of  measurements 


3. Model  identification 
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The  ultimate  goal  of  any  analysis  technique  is  the  complete  characterization  of  the 
system  under  investigation.  In  this  thesis  it  is  proposed  that,  the  basic  building  blocks  on  which 
any  such  effort  is  based  are  experimental  protocol,  measurement  model,  fit  criteria  and  model 
identification.  Examination  of  the  flow  diagram  in  Figure  2 shows  that  the  there  is  a very  strong 
interaction  between  all  the  key  components.  The  whole  process  is  iterative,  where  the  starting 
point  is  the  prior  knowledge  of  the  system. 

1.1.1  Optimal  Design  of  Experiments 

A thorough  understanding  of  the  theory  of  impedance  spectroscopy  and  the  working 
of  the  instrumentation  is  crucial  for  the  optimal  design  of  impedance  experiments.  In  chapter  2, 
the  principles  of  impedance  measurements  are  discussed.  A brief  review  of  the  instruments 
used  to  collect  impedance  data  is  presented.  The  software  to  collect  impedance  data  was 
developed  in-house.  The  details  of  the  software  are  also  presented  in  chapter  2.  Some  of  the 
unique  features  of  the  software  are  the  automatic  selection  of  measuring  resistor  during  the 
course  of  the  experiment  and  the  algorithm  for  galvanostatic  measurements.  The  importance  of 
optimal  choice  of  measuring  resistor  was  realized  when  the  analysis  of  the  data  by 
measurement  model  indicated  that  improper  choice  of  the  measuring  resistor  makes  a 
significant  contribution  to  the  noise  in  the  measurement. 

In  addition  to  the  proper  design  of  software  for  to  run  the  experiments,  there  are  other 
issues  which  arise  due  to  the  inability  of  impedance  spectroscopy  to  serve  as  a stand-alone 
method  for  identification  of  a correct  model.  These  issues  can  be  addressed  experimentally  by 
including  additional  analysis  techniques  or  by  incorporating  multiple  or  more  directed  forcing 
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functions.  Discussion  of  other  analysis  techniques  and  the  theory  behind  the  optimal  design  of 
experiments  are  outside  the  scope  of  this  thesis. 

1.1.2  Characterization  of  Measurements 

The  issue  of  identification  of  the  characteristics  of  impedance  data  is  particularly 
important  when  model  selection  is  guided,  at  least  in  part,  by  regression  of  the  model  to 
experimental  data.  In  this  work,  a “measurement”  model  is  used  as  a statistical  tool  to  identify 
fundamental  properties  of  impedance  measurements.  A measurement  model  is  defined  as  being 
a superposition  of  a series  of  lineshapes,  and,  with  a sufficient  number  of  terms,  it  provides  a 
statistically  adequate  fit  to  any  consistent  data  set.  The  concept  of  the  measurement  model  is 
discussed  in  chapter  3.  It  is  shown  that  a measurement  model  based  on  Voigt  circuit  elements 
can  adequately  represent  most  electrochemical  systems.  FORTRAN  code  for  building  a 
measurement  model  based  on  Voigt  circuit  elements  is  shown  in  Appendix  A 

Selection  of  an  appropriate  model  for  a given  system  requires  an  ability  to  identify  the 
portion  of  a spectrum  corrupted  by  instrumental  anomalies  or  by  non-stationary  behavior.  The 
Kramers-Kronig  relations  [5]  provide  a unique  transformation  that  can  be  used  to  predict  one 
component  of  the  impedance  if  the  other  is  known  over  the  frequency  limits  of  zero  to  infinity 
and  if  the  data  satisfy  the  constraints  that  the  system  was  causal,  linear,  stable,  and  (by 
implication)  stationary.  The  Kramers-Kronig  relations  can,  in  principle,  be  applied  with  a 
suitable  extrapolation  of  the  data  into  the  unmeasured  frequency  domain.  While  extrapolation 
approaches  have  been  applied  to  some  experimental  data  with  success,  any  approach  toward 
extrapolation  can  be  applied  over  only  a small  frequency  range  and  cannot  be  applied  at  all  if 
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the  experimental  frequency  range  is  so  small  that  the  data  do  not  show  a maximum  in  the 
imaginary  impedance. 

The  limitations  associated  with  extrapolation  of  impedance  spectra  are  resolved  in  this 
work  by  taking  advantage  of  the  lack  of  uniqueness  of  impedance  models.  Since  the  model  has 
been  shown  to  provide  an  adequate  representation  of  a consistent  spectrum,  failure  to  fit  the 
data  can  be  attributed  to  inconsistency  with  the  Kramers-Kronig  relations.  In  effect  the  problem 
of  extrapolation  is  converted  into  a problem  of  regression  of  the  measurement  model  to  the 
experimental  data.  The  issues  concerning  consistency  of  data  with  Kramers-Kronig  relations 
and  an  algorithm  based  on  Monte-Carlo  simulation  are  discussed  in  chapter  5.  The  FORTRAN 
code  for  the  algorithm  is  shown  in  Appendix  B. 

From  the  previous  sections,  it  is  clear  that  regression  plays  a very  crucial  role  in 
building  of  the  measurement  model  for  the  system,  in  the  algorithm  to  determine  consistency  of 
impedance  data  and  finally  in  identifying  a process  model.  It  should  be  noted,  however,  that  the 
impedance  data  exhibit  strong  heteroscedasticity  (non-uniform  variance)  [8]  and  the  models  are 
usually  highly  non-linear.  Hence,  regression  is  sensitive  to  the  weighting  applied  to  the  data. 
For  non-linear  regression  with  heteroscedasticity,  weighting  with  the  variance  of  the  data  is 
preferred.  In  chapter  4 the  error  structure  for  impedance  data  is  discussed,  and  a model  for  the 
variance  of  impedance  data  is  proposed.  Identification  of  the  error  structure  is  critical,  not  only 
because  the  variance  of  the  data  is  the  preferred  weighting  for  regression,  but  because  the  error 
in  the  experiment  indicates  a limit  as  to  how  well  a model  can  fit  data  and  because  knowledge 
of  the  error  structure  can  guide  improvement  in  the  experimental  design. 
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1.1.3  Model  Identification 

As  stated  earlier,  the  ultimate  aim  of  an  analysis  technique  is  the  characterization  of  the 
system  under  investigation  by  a process  model.  A process  model  is  a deterministic  model  built 
upon  hypothesized  physical  and  chemical  phenomena.  In  chapter  6 the  approach  outlined  in 
Figure  2 is  applied  to  a metal  hydride  system.  A process  model  for  the  system  is  developed  and 
the  inadequacies  of  the  model  are  shown. 

In  chapter  7,  the  approach  is  applied  to  electro-hydrodynamical  impedance  of  a 
rotating  platinum  disc  electrode.  It  is  shown  that  the  approach  proposed  in  this  thesis  helps  in 
the  characterization  of  the  system  by  improving  the  regression  of  the  process  model  to 
experimental  data. 


1.2  New  Directions  for  Statistical  Analysis 

The  measurement  model  appears  to  provide  a useful  framework  for  a statistical  analysis 
of  impedance  data.  The  approach  is  currently  limited  to  evaluating  consistency  of  data  with  the 
Kramers-Kronig  relations,  to  filtering  data  for  lack  of  replicacy  (a  step  in  the  experimental 
identification  of  the  error  structure  of  impedance  measurements),  and  to  providing  a rough 
guide  for  model  identification.  Work  is  needed  to  improve  the  use  of  measurement  models  for 
model  identification,  for  characterization  of  industrial  data  (where  extrapolation  is  often  more 
important  than  model  identification,  and  for  characterization  of  non-stationary  phenomena. 
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Figure  1.1:  Flow  diagram  for  the  measurement  and  characterization  of  a material-electrode 
system  [7]  (reproduced  with  permission  from  the  author). 


Figure  1.2:  Flow  diagram  for  impedance  spectroscopy  proposed  in  this  work. 


CHAPTER  2 


IMPEDANCE  SPECTROSCOPY:  EXPERIMENTAL  DESIGN 

In  this  chapter,  the  design  of  an  impedance  experiment  is  outlined.  The  theory  of 
impedance  measurements  using  a single-sine  technique  is  summarized.  The 
instrumentation  required  to  run  an  impedance  experiment  is  discussed.  A software  to  run 
the  impedance  was  developed  in  house  using  Labview®  programming  language.  The  key 
features  of  the  software,  which  differentiate  it  from  the  commercially  available  software, 
are  highlighted. 

2 . 1 Theory  of  Impedance  Measurement 
The  fundamental  approach  of  all  impedance  methods  is  to  apply  a small  amplitude 
sinusoidal  excitation  signal  (voltage  or  current)  to  the  system  under  investigation  and 
measure  the  response  (current  or  voltage)  [4,7,10,11],  The  excitation  signal  can  be  applied 
in  several  ways.  The  two  most  commonly  used  methods  are  the  multi-sine  technique  [12, 
13]  and  the  single-sine  technique.  Fast  measurement  time  and  mild  perturbation  of  the 
system  under  investigation  have  been  stated  to  be  the  strength  of  the  multi-sine  technique, 
though  it  has  been  shown  by  other  authors  [4]  that  using  a fast  frequency  response 
analyzer  the  difference  in  the  measurement  time  for  a single-sine  and  a multi-sine 
technique  is  small.  One  of  the  disadvantages  of  the  technique  is  that  it  has  a small 
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frequency  range.  The  multi-sine  technique  is  also  sensitive  to  harmonic  distortion  and  is 
vulnerable  to  aliasing. 

In  this  work  impedance  data  was  collected  using  the  single-sine  technique.  The 
theory  underlying  a single-sine  impedance  measurement  will  be  discussed  in  detail 
subsequent  sections. 


2.1.1  Single-Sine  Technique 


The  single-sine  technique  uses  a digital  correlation  frequency  response  analyzer  to 
generate  the  excitation  waveform  (usually  a sine  wave  with  a single  frequency)  and 
analyze  the  response  to  determine  the  impedance  of  the  system.  In  Figure  2. 1 a non-linear 
I-V  curve  for  a theoretical  electrochemical  system  is  shown  [4],  A low  amplitude  sine 
wave  A£sin(cot),  of  a particular  frequency,  is  superimposed  on  the  dc  polarization  voltage 
E0.  This  results  in  a current  response  of  a sine  wave  A/sin(cot+4>)  superimposed  on  the  dc 
current  I0.  The  Taylor  series  expansion  for  the  current  is  given  by 
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If  the  magnitude  of  the  perturbing  signal  A E is  small,  then  the  higher  order  terms 
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AE2+....  in  equation  2.1  can  be  assumed  to  be  negligible.  The  impedance  of 
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the  system  can  then  be  calculated  as  . This  calculation  is  done  at  several  frequencies  to 


get  the  impedance  of  the  system  as  a function  of  frequency.  In  general,  the  response 
includes  the  desired  information  at  the  required  frequency,  corrupted  by  harmonics  and 


noise. 
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A potentiostat  is  used  to  measure  and  control  the  potential  and  current.  The 
sinusoidal  perturbation  is  applied  by  a Frequency  Response  Analyzer  (FRA).  The  principle 
of  the  measurement  of  the  transfer  function  using  a FRA  is  schematized  in  Figure  2.2 
[4,14],  Consider  a sinusoidal  perturbing  voltage  signal  V(t),  represented  as 

V(t)  = V0  sin(rttf)  (2.2) 

where  Vo  is  the  amplitude  and  co  is  the  frequency.  The  current  response  1(f)  can  be 
represented  as 

I(t)  - sin[°fr  + <!>(©)]  + 2 sin(nmt-§m)+n(t)  (2.3) 

m 

where  ^T(co)ei'K“)  is  the  cell  transfer  function  and  <J>(co)  is  the  phase  shift  between  the 
perturbing  function  and  the  cell  response.  The  first  term  on  the  right  side  of  equation  2.3  is 
the  fundamental  component  of  the  cell  response,  the  second  term  is  the  sum  of  the  various 
higher  order  harmonics  and  the  third  term  is  the  random  noise.  As  shown  in  Figure  2.2  the 
cell  response  I(t)  is  correlated  with  two  synchronous  reference  signals,  one  in  phase  with 
V(t)  and  the  other  90°  out  of  phase,  i.e.  sin(cot)  and  cos(cot),  in  order  to  calculate  the  real 
and  imaginary  components  of  the  transfer  function,  as 


I'(g>)  = J/(0sin(cfr)^ 

/"(to)  = ~ J/(f)  cos(cof)fifr 
1 0 


(2.4) 


where,  / and  / are  the  real  and  imaginary  part  the  response  in  the  frequency  domain. 
Substitution  of  equation  2.2  in  equation  in  equation  2.3  results  in 
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and 


r'(co)  = V0K(a)j  sin[otf  + <|>(cd)]  sin(co/)^ 

0 

+ 7 1 Z sin[/M0f  - (f)  m ] sin(afr)  dt 

* 0 m 

+ip]n(t)sin((Qt)dt 

* o 

r'(co)  = F0  A'(co) } sin[cot  + 4>(ffl)]  cos(ofr)  dt 

0 

sin[/w®f  — 4>m ] cos(cof) dt 

* 0 m 

+-7z]n(t)  cos((o  t)dt 


(2.5) 


(2.6) 


If  the  noise  n(t)  is  completely  random  (uncorrelated),  then,  as  the  integration  time  T goes 
to  infinity,  the  last  integrals  in  equations  2.4  and  2.5  go  to  zero.  Another  way  of  describing 
this  result  is  that  the  FRA  acts  as  a filter  with  an  infinitely  narrow  bandwidth,  centered 
around  the  generator  frequency.  However  the  integration  time  T cannot,  in  practice  be 
infinite;  hence  the  actual  bandwidth  of  the  filter  is  dependent  on  T.  If  the  integration  is 
carried  out  over  Nf  periods  of  sinusoidal  perturbation,  the  equivalent  filter  selectivity  is 
given  by 

A/  1 

(2-7) 

where  fi  is  the  center  or  the  analyzed  frequency  and  A f is  the  bandwidth  [4],  The 
contribution  of  the  higher  harmonic  terms,  represented  by  the  second  term  in  equations  2.4 
and  2.5,  go  to  zero  if  the  integrals  are  carried  over  multiples  of  2n  due  to  the  orthogonal 
properties  of  sine  and  cosine  functions.  Therefore,  the  integral  in  equations  2.4  and  2.5 


can  be  evaluated  as 
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/'(to)  = — J I(t)  sin(co t)dt  = V0K((£>)  cos(|)(q))  as  T ->  oo; 

0 

r (2.8) 

/'(CO)  = — J I(t)cos((Qt)dt  = F0/f(co)sin<j>(co)  as  T — > oo. 

0 

A two  channel  FRA  measures  the  in-phase  and  out-of-phase  components  of  a signal  with 
respect  to  the  perturbing  signal.  If  the  measurement  is  done  in  potentiostatic  mode,  then 
the  measured  current  is  given  by 

I {a)  = Gj  RKi  (VqK(co)  cos<p(co)  + jV0K(a>)  sin  </>(a>))  (2.9) 

and  the  measured  potential  is  given  by 

V(co)  = GvRKlV0  (2.10) 

where  K\  is  the  transfer  function  of  the  potentiostat,  R is  a standard  resistor  for  measuring 
current,  and  Gv  and  G\  are  the  gains  of  the  voltage  and  current  amplifiers.  The 
potentiostats  are  designed  such  that  Gv-Gr=l  at  low  frequencies  and  they  are  balanced  at 
high  frequencies,  i.e.  they  have  the  same  phase  shifts  at  a given  frequency.  Therefore, 


Z(q)  = 


n«p>_ gvrkxv0 

/(to)  G,  R Kx  (yoK(<x>)  cos<j)(co)  + jK{ co)  sin  <{>(©)) 

1 


(2.11) 


K(a)e 


>«K“) 


where  Z(co)  is  the  impedance  of  the  cell. 


2.2  Instrumentation 

A typical  single-sine  experimental  setup  includes  an  electrochemical  cell,  which  is 
the  system  under  investigation.  The  potential  is  controlled  by  a potentiostat.  A frequency 
response  analyzer  applies  the  sine  wave  and  analyzes  the  response  of  the  system  to 
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determine  the  impedance  of  the  system.  The  instrumentation  is  interface  to  a high  speed 
digital  computer,  which  controls  the  experiment  and  stores  the  data. 

2.2.1  The  Electrochemical  Cell 

The  electrochemical  cell  in  an  impedance  experiment  can  consist  of  two,  three,  or 
four  electrodes.  The  most  basic  form  of  the  cell,  shown  in  Figure  2.3a,  has  two  electrodes 
[11].  Usually  the  electrode  under  investigation  is  called  the  working  electrode,  and  the 
electrode  necessary  to  close  the  electrical  circuit  is  called  the  counter  electrode.  The  cell 
can  be  include  an  electrolyte  solution  or,  for  solid-state  system,  may  not  have  any 
electrolyte.  A two-electrode  configuration  for  the  cell  is  used  when  the  precise  control  of 
the  potential  is  not  very  critical.  This  arrangement  is  used  to  investigate  electrolyte 
properties,  such  as  conductivity,  or  to  characterize  solid  state  systems. 

The  three-electrode  configuration  for  an  electrochemical  cell  is  most  common  for 
typical  electrochemical  applications.  A schematic  of  such  a cell  is  shown  in  Figure  2.3b.  A 
third  electrode  (the  reference  electrode)  is  used  to  determine  or  control  the  potential  of  the 
working  electrode  precisely.  Since  the  absolute  potential  of  a single  electrode  cannot  be 
measured,  all  potential  measurements  in  an  electrochemical  systems  are  performed  with 
respect  to  a reference  electrode.  A reference  electrode,  therefore,  should  be  reversible,  and 
its  potential  should  remain  constant  during  the  course  of  the  measurement. 

Finally,  a four-electrode  cell,  shown  in  Figure  2.3c,  is  used  to  analyze  processes 
occurring  within  the  electrolyte,  between  two  measuring  electrodes  separated  by  a 
membrane.  In  this  configuration,  the  working  electrode  and  the  counter  enable  current 
flow.  This  kind  of  a cell  is  usually  used  to  study  ion  transport  through  a membrane. 
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2.2.2  Potentiostat 

Electrochemical  systems  require  a precise  control  of  the  potential  of  the  working 
electrode  with  respect  to  the  reference  electrode  because  the  rates  of  electrochemical 
reactions  depend  on  the  potential  of  the  electrodes.  A potentiostat  is  therefore  a very 
important  instrument  for  any  electrochemical  impedance  measurement.  When  evaluating  a 
potentiostat,  one  looks  at  three  main  components,  the  power  amplifier,  the  electrometer 
and  the  I/E  converter  [12,15], 

The  power  amplifier  of  a potentiostat  controls  the  counter  to  working  electrode 
voltage  to  achieve  the  proper  reference  to  working  electrode  voltage.  Any  limitations  in 
the  frequency  response  of  the  power  amplifier  will  adversely  affect  the  control  of  desired 
voltage  between  the  working  electrode  and  the  reference  electrode.  The  electrometer  in  a 
potentiostat  is  used  to  sense  the  potential  difference  between  the  working  electrode  and 
the  reference  electrode  and  feed  it  back  to  the  power  amplifier.  This  feedback  control  of 
the  potential  can  result  in  oscillations  and  can  thus  influence  the  overall  stability  of  the 
electrochemical  system.  Different  cells  with  dissimilar  electrical  characteristics  have 
varying  effects  on  the  stability  of  the  system.  To  maintain  an  oscillation  free  system,  the 
cell  differences  must  be  compensated  for  by  changes  in  the  rest  of  the  system. 
Potentiostats  allow  a choice  of  bandwidth  appropriate  for  different  cell  configurations. 
Selection  of  a suitable  bandwidth  limits  the  range  of  frequencies  at  which  these  oscillations 
can  occur  and  can  hence  help  stabilize  the  system.  Finally,  the  measurement  of  current  in 
an  experiment  is  done  by  passing  the  current  through  a known  resistor  and  measuring  the 
voltage  drop  across  it.  A potentiostat  typically  has  a wide  range  of  measuring  to  choose 
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from.  The  selection  of  the  measuring  resistor  depends  on  the  magnitude  of  the  current 
being  measured. 

Two  different  potentiostats  were  used  in  this  work.  The  solid  state  data[16]  and 
copper  corrosion  data  [17]  were  collected  with  the  SI 1286  potentiostat.  The  current 
measuring  resistors  available  on  the  SI1286  potentiostat  are  tabulated  in  Table  2.1  [18]. 
The  SI1286  offers  up  to  ten  bandwidths  (A-J)  to  ensure  the  stability  of  the  system  under 
investigation  [18],  A PAR  273 A potentiostat  was  used  for  measuring  the  impedance  of 
metal  hydrides.  The  measuring  resistors  offered  on  the  PAR  273A  potentiostat  are  shown 
in  Table  2.2  [19],  The  PAR  273A  potentiostat  provides  the  choice  of  two  bandwidths 
(high  speed  and  high  stability). 

2.2.3  Frequency  Response  Analyzer 

The  Solartron  1250  Frequency  Response  Analyzer  was  used  for  all  the  impedance 
experiments  presented  in  this  work.  The  function  of  the  FRA  can  be  divided  into  two 
major  components,  generator  and  analyzer.  The  FRA  is  capable  of  generating  sine,  square 
or  triangle  waveform  at  frequencies  from  lOpHz  to  65,500Hz.  These  waveforms  are 
synthesized  with  a 10MHz  crystal  oscillator.  The  analyzer  section  of  the  FRA  can 
simultaneously  measure  and  analyze  two  input  channels  [20], 

2.2.4  Computer  Interface 

The  potentiostat  and  the  frequency  response  analyzer  were  interfaced  to  a high 
speed  personal  computer  through  a IEEE  488  GPIB  interface.  The  impedance  experiments 
were  performed  using  the  ZPLOT®  software  [21]  and  a software  developed  in-house, 
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using  the  Labview®  programming  language.  The  details  of  the  in-house  software  are  given 
in  the  next  section. 

2.3  Software 

The  impedance  experiments  conducted  with  the  SI1286  potentiostat  and  the  1250 
FRA  were  conducted  using  the  ZPLOT®  software.  In  the  course  of  the  work  it  was 
discovered  that  the  optimum  choice  of  the  measuring  resistor  is  crucial  to  minimize  the 
errors  in  the  experimental  data.  The  ZPLOT®  software  allows  the  use  of  only  two 
measuring  resistors  during  the  course  of  an  impedance  scan.  These  measuring  resistors 
have  to  be  fixed  before  the  start  of  the  experiment.  A need  to  have  a software  which  can 
automatically  select  an  optimum  measuring  resistor  was  felt.  None  of  the  commercially 
available  software  had  this  feature.  To  overcome  the  problem  of  limiting  the  experimental 
parameters  to  the  available  options  on  commercial  software  an  in-house  software  to  run 
impedance  experiments  was  developed.  The  Labview®  language  was  chosen  to  write  the 
software.  The  Labview®  software  provides  a user  friendly  interface  and  was  found  to  be 
very  powerful  for  the  application. 

The  software  was  developed  for  impedance  experiments  to  be  conducted  with  a 
PAR273A  potentiostat  in  conjunction  with  a Solartron  1250  FRA.  The  software  provides 
a user  friendly  menu  driven  interface. 


18 


2.3.1  Menu  Options 

In  this  section  the  various  menu  options  are  discussed  in  detail.  Proper 
understanding  of  the  potentiostat  and  FRA  settings  is  crucial  to  optimize  the  experimental 
protocol. 

2 . 3 . 1 . 1 Potentiostat  Settings 

Pot.  or  Galv.  Mode:  The  potentiostat  can  be  operated  in  the  potentiostatic  or 
galvanostatic  mode.  In  the  potentiostatic  mode,  the  user  specifies  the  potential  bias  (see 
Sweep  Bias  menu)  to  be  applied  to  the  cell.  The  potential  of  the  working  electrode  with 
respect  to  the  reference  electrode  is  controlled  by  the  potentiostat.  A potential 
perturbation  is  applied  to  the  cell,  and  the  resulting  current  is  measured  to  determine  the 
impedance  of  the  system. 

In  the  galvanostatic  mode,  the  user  specifies  the  bias  current  to  be  maintained 
during  the  course  of  the  experiment  (see  Sweep  Bias  menu).  The  potentiostat  maintains 
the  desired  current  by  controlling  the  potential  difference  between  the  counter  electrode 
and  the  working  electrode. 

Bandwidth:  The  PAR  273A  potentiostat  provides  a choice  of  two  bandwidths, 
high  stability  and  high  speed.  In  the  high  speed  mode,  the  potentiostat  provides  the  fastest 
response.  This  setting  is  recommended  as  the  default,  as  it  provides  the  widest  frequency 
range.  This  setting  can  cause  ringing  and  oscillations  with  cells  that  have  large  capacitance 
and  low  resistance.  The  high  stability  mode  provides  extreme  stability  at  somewhat 
reduced  speed.  The  choice  of  this  setting  reduces  the  frequency  range  of  the  impedance 
experiment. 
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Current  Range:  The  eight  current  measuring  resistors  available  on  PAR273 A are 
listed  in  table  5.2.  In  the  commercially  available  software,  the  user  chooses  the  measuring 
resistors  to  be  used  during  the  experiment  at  the  start  of  the  experiment.  Usually  the  user 
is  allowed  a choice  of  one  or  two  measuring  resistor  for  the  entire  scan.  For  systems 
where  the  impedance  changes  by  several  order  of  magnitude  over  the  course  of  the 
experiment,  this  could  result  in  a low  signal  to  noise  ratio.  To  overcome  this,  one  usually 
breaks  the  frequency  range  into  several  parts;  conducts  separate  experiments  and 
concatenates  the  data  later.  This  approach  could  lead  in  an  increase  in  the  time  required  to 
collect  data.  Also,  the  choice  of  the  measuring  resistor  is  arbitrary.  In  this  work,  the 
program  checks  to  see  if  the  proper  measuring  resistor  is  being  used  and  then  changes  it 
automatically  when  appropriate. 

There  are  several  probable  strategies  to  select  an  optimum  measuring  resistor.  The 
absolute  value  of  the  impedance  at  the  frequency  at  which  the  data  is  to  be  collected  is  the 
crucial  parameter  for  any  algorithm  to  select  a measuring  resistor  as  it  determines  the  total 
current.  All  strategies  would  have  to  use  some  kind  of  predictive  algorithm  to  guess  the 
absolute  value  of  the  impedance  for  the  next  measurement  based  on  previous 
measurements.  One  possible  strategy  would  be  to  do  a polynomial  fit  to  the  data  that  have 
already  been  collected  and  to  extrapolate  to  predict  the  next  value.  The  next  value  can  also 
be  predicted  by  the  regression  of  a measurement  model  (to  be  discussed  in  the  next 
chapter)  to  the  data  and  extrapolation  of  the  measurement  model.  This  would  provide  a 
more  robust  prediction.  The  problem  with  these  strategies  is  that  they  would  increase  the 
data  collection  time.  The  algorithm  described  here  uses  the  current  value  for  the 
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impedance  as  the  guess  for  the  next  data  point.  The  algorithm  described  here  is  valid  for 
impedance  experiments  where  the  absolute  value  of  the  impedance  increases  with 
frequency  (this  is  usually  true  when  the  frequency  sweep  is  from  high  frequency  to  low 
frequency).  Although  this  algorithm  can,  in  some  cases,  result  in  inaccurate  prediction  (for 
cases  where  there  is  very  large  change  in  the  impedance  of  the  system  from  one  frequency 
to  the  next),  it  was  found  to  work  well  for  the  systems  considered  in  this  work.  The 
algorithm  did  not  result  in  any  significant  increase  in  the  data  collection  time.  The 
algorithm  is  as  follows: 

1.  At  the  beginning  of  the  experiment  the  user  specifies  a current  range  (measuring 
resistor).  Prior  knowledge  of  the  system  is  used  to  provide  an  appropriate  current 
range. 

2.  Measurement  at  the  first  frequency  of  the  impedance  scan  is  done  using  the 
measuring  resistor  provided  by  the  user.  The  absolute  value  of  the  resultant 
impedance  is  stored. 

3.  The  impedance  value  is  used  to  calculate  the  expected  current  for  the  next 
measurement. 

4.  The  expected  current  is  compared  with  the  current  range  (measuring  resistor). 
If  the  calculated  current  is  smaller  than  the  appropriate  current  range  then  a 
command  is  sent  to  the  potentiostat  to  switch  to  the  proper  current  range.  For 
example,  if  the  current  range  being  used  is  1mA  to  10mA  and  the  measured 
current  is  .9  mA  then  the  measuring  resistor  is  changed  such  that  the  current  range 


is  ,1mA  to  1mA. 
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5.  The  smallest  available  current  range  for  PAR  273 A is  InA-lOOnA.  If  the 
expected  current  is  lower  that  InA,  then  the  algorithm  fixes  the  measuring  resistor 
to  the  largest  available  value  of  10MQ. 

It  is  important  to  choose  a proper  measuring  resistor  at  the  beginning  of  the  experiment 
such  that  there  is  no  current  overload;  otherwise,  the  algorithm  will  fail. 

2.3. 1.2  Sweep  Bias 

Run  Number  (Bias):  The  software  allows  the  user  to  run  a multiple  set  of  scans. 
This  setting  is  used  to  set  the  bias  for  each  individual  scan.  For  experiments  under 
potentiostatic  control,  the  user  specifies  the  bias  potential  to  applied  for  each  run.  This 
bias  or  polarization  voltage  can  be  specified  vs.  the  reference  electrode  (this  is  the  default) 
or  vs.  the  open  circuit  potential.  As  the  potentiostat  applies  a potential  with  respect  to  a 
reference  electrode,  if  the  user  chooses  the  “vs.  open  circuit  potential”  option,  then  the 
potential  specified  is  added  to  the  open  circuit  potential,  measured  at  the  beginning  of  each 
scan.  This  way  the  software  ensures  a correct  application  of  polarization  voltage.  For 
impedance  experiment  under  galvanostatic  control,  the  bias  current  is  specified  in  mA. 

Sweep  Run  Mode:  The  frequency  sweep  can  be  applied  in  four  different  ways.  In 
the  log-up  mode  the  frequencies  are  spaced  logarithmically,  and  the  sweep  starts  at  the 
low  frequency  end  and  ends  at  the  high  frequency  end.  In  the  log-down  mode,  the  sweep 
starts  at  the  high  frequency  end  and  ends  at  the  low  frequency  end.  The  software  uses  this 
as  the  default  mode.  In  the  linear-up  and  linear-down  modes  the  frequencies  are  spaced 


linearly. 
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2.3 . 1 .3  Frequency  Sweep  Parameters 

Maximum  Frequency:  The  user  specifies  the  maximum  frequency  (in  kHz)  for  the 
frequency  scan.  The  maximum  allowable  frequency  for  the  FRA  1250  is  65.535kHz. 

Minimum  Frequency:  The  user  specifies  the  minimum  frequency  (in  Hz)  for  the 
frequency  scan.  The  minimum  allowable  frequency  for  the  FRA  1250  is  1.0  x 10'5Hz. 

Sweep  Increment:  If  log  increment  is  chosen  then  the  user  specifies  the  number  of 
steps  per  decade  for  the  frequency  scan.  If  linear  increment  is  chosen  then  the  user 
specifies  the  number  of  steps  into  which  the  whole  frequency  range  is  to  divided. 

2.3. 1.4  FRA  Display 

Source:  The  FRA  has  two  available  channels.  For  electrochemical  impedance 
experiments,  one  channel  is  used  to  measure  the  voltage  signal  and  the  other  to  measure 
the  current  signal.  The  source  setting  controls  the  display  on  the  front  panel  of  the  FRA. 
The  choices  available  are  Chi,  Ch2,  Chl/Ch2  and  Ch2/Chl.  In  this  work  channel  1 is 
used  to  measure  the  current  and  channel  2 is  used  to  measure  the  voltage.  Hence  the 
default  setting  used  is  Ch2/Ch2  which  is  the  impedance  of  the  cell. 

Co-ordinates:  The  FRA  allows  the  displays  to  be  in  Cartesian,  polar  or  log  polar 
co-ordinates.  The  default  used  here  is  the  Cartesian  co-ordinate.  The  polar  co-ordinates 
are  calculated  in  the  program  and  are  plotted  along  with  the  plots  in  the  Cartesian  co- 
ordinate system. 

Error:  The  FRA  produces  an  error  code  and  a beep  if  the  instrument  detects  any 
kind  of  error  (The  listing  of  the  error  codes  is  given  in  the  FRA  manual).  This  error  beep 
can  be  turned  off  if  the  user  so  desires.  The  default  setting  is  on. 
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2.3. 1.5  FRA  Settings 

Integration  Time:  The  user  can  choose  the  integration  time  in  seconds  or  cycles. 
The  minimum  values  allowed  are  10ms  or  1 cycle(  whichever  is  longer)  of  the  fundamental 
frequency  for  frequencies  up  to  655Hz.  For  frequencies  greater  than  655Hz,  10ms  or  61 
cycles  (whichever  is  longer)  is  used.  If  the  integration  time  is  chosen  in  seconds  then  the 
FRA  automatically  rounds  off  the  number  to  give  the  exact  number  of  cycles  at  the  given 
frequency.  The  integration  time  is  chosen  such  that  there  is  a compromise  between  the 
speed  of  measurement  and  acceptability  of  small  errors  in  the  measurement.  For 
measurements  where  the  random  noise  level  is  expected  to  high,  the  number  of  cycles  or 
time  needs  to  be  increased  to  improve  the  signal  to  noise  ratio.  If  the  integration  time  is 
too  small,  then  the  FRA  gives  error  82. 

Auto-Integration:  Usually  it  is  difficult  for  the  user  to  decide  how  long  the 
integration  time  should  be.  In  such  a case,  the  user  sets  a large  integration  time  and  lets 
the  FRA  use  the  auto-integration  mode.  In  the  auto-integration  mode  the  measurement  is 
continued  until  either  the  variation  in  the  running  average  becomes  sufficiently  small  or  the 
maximum  time  allowed  has  expired.  There  are  two  levels  of  acceptable  variation  available 
on  the  FRA  1250,  short  or  long.  The  reads  are  averaged  until  the  standard  deviation  falls 
below  the  following  limits: 

short:  = ± 10%  of  reading  ± 0.01%  of  frill  scale, 
long:  = ± 1.0%  of  reading  ± 0.001%  of  full  scale. 

Student ‘t’  test  is  applied  to  the  standard  deviation  to  ensure  a 90%  confidence  level. 
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Measurement  Delay:  This  allows  a time  delay  to  be  interposed  between  applying 
the  generator  waveform  to  the  cell  and  taking  a reading.  This  allows  the  cell  to  settle  after 
the  input  frequency  is  changed  during  a sweep.  The  delay  can  be  entered  as  time  in 
seconds  or  as  a specified  number  of  cycles.  The  FRA  converts  the  delay  setting  into  time 
value  in  seconds  before  applying  it  to  the  system  under  test. 

Harmonic:  The  FRA  can  be  used  to  analyze  the  response  of  the  system  under 
investigation  at  any  harmonic  of  the  fundamental  up  to  the  16th.  For  single  sine 
experiment  only  the  fundamental  harmonic  is  used. 

Range:  The  FRA  allows  five  selections  of  voltage  ranges  for  the  two  input 
channels  (30  mV,  300  mV,  3 V,  30  V,  300  V)  and  an  Auto  range  selection.  If  the  values 
of  the  voltages  expected  is  known  then  the  proper  range  should  be  used  to  get  the  best 
resolution,  otherwise  Auto  range  should  be  used. 

Coupling:  There  are  two  modes  of  coupling  the  signal  to  the  analyzer  input,  dc 
and  ac.  The  dc  mode  is  recommended  as  it  introduces  the  minimum  phase  shift.  If  the 
signal  is  suspected  to  have  been  corrupted  by  unwanted  dc  component  then  ac  coupling 
should  be  used  to  block  the  ac  component. 

Input:  The  input  to  the  FRA  can  be  routed  via  the  4mm  sockets  on  either  the  front  or  the 
rear  panel.  The  user  has  to  specify  the  location  of  the  input. 

2.3. 1.6  FRA  File  Settings 

Filename:  The  user  specifies  the  name  of  the  file  into  which  data  is  to  written.  The 
filename  should  not  be  longer  than  six  characters.  The  program  automatically  appends 
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rdat  to  the  file  in  which  the  impedance  data  is  stored,  where  / stands  for  the  run  number. 
The  open  circuit  information  is  stored  in  filename,  oc. 

Directory:  The  user  can  specify  the  directory  in  which  the  file  is  to  be  written. 

Overwrite:  If  this  is  set  as  yes  then  the  old  file  is  overwritten. 

2.3. 1.7  Generator  Setting 

Waveform : The  FRA  allows  a choice  of  three  waveforms:  sine,  triangle,  and 
square.  The  waveform  used  in  this  work  is  the  sine  wave. 

Stop  at:  This  setting  enables  the  generator  to  be  stopped  at  0°,  90°,  180°,  270° 
point  of  the  waveform.  This  facility  ensures  that  the  waveform  always  start  from  a known 
point.  The  default  setting  of  0°  enables  any  offset  in  the  system  under  test  to  be  adjusted 
to  zero. 

Scan  Mode:  The  two  available  modes  are  sweep  and  single  frequency.  Under  the 
sweep  mode  impedance  is  measured  for  the  full  frequency  range.  If  single  frequency  mode 
is  chosen  then  the  measurement  is  done  only  at  one  frequency. 

Amplitude:  For  potentiostatic  mode  the  user  specifies  the  perturbation  potential 
amplitude  in  mV.  The  amplitude  of  perturbation  should  small  enough  to  ensure  linearity 
and  large  enough  to  get  a good  signal  to  noise  ratio. 

For  performing  impedance  experiments  in  the  galvanostatic  mode,  the  traditional 
approach  has  been  to  choose  a fixed  perturbation  current  as  it  is  done  for  the 
potentiostatic  case.  The  problem  with  this  approach  is  that  typically  the  impedance  of  the 
system  under  investigation  increases  with  the  decrease  in  frequency.  This  results  in  the 
increase  in  the  voltage  required  to  maintain  a constant  current  perturbation,  which  could 
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lead  to  the  system  being  non-linear.  To  overcome  this  problem  of  non-linearity  a new 
approach  is  used  in  this  work.  This  approach  takes  advantage  of  the  fact  that  for  linear 
system  the  response  is  independent  of  the  magnitude  of  perturbation.  Instead  of  applying  a 
constant  current  perturbation,  it  is  varied  in  such  a way  that  the  voltage  required  to  drive 
the  current  is  always  kept  below  a target  voltage.  The  user  specifies  a target  voltage 
amplitude  to  ensure  linearity.  At  each  frequency  the  absolute  value  of  the  impedance  of  the 
system  is  calculated.  The  ratio  of  the  target  voltage  amplitude  to  the  calculated  impedance 
determines  the  value  of  the  current  perturbation  to  be  applied  at  the  next  frequency,  i.e., 
v 

target  . . 

|z((o)|  — (®)  perturbation  (2-12) 

where  5.12  V target  is  the  user  specified  target  amplitude,  I Z(co)|  is  the  absolute  value  of  the 
impedance  at  a certain  frequency  and  I(a> )perturbation  is  the  calculated  perturbation  for  the 
next  frequency.  At  the  beginning  of  the  experiment  the  user  specifies  a guessed  impedance 
value  to  calculate  the  current  perturbation  at  the  first  frequency. 

2.3.2  Typical  Experimental  Results 

The  front  panel  for  running  an  impedance  experiment  in  the  potentiostatic  mode  is 
shown  in  Figure  2.5,  and  a front  panel  for  the  galvanostatic  mode  is  shown  in  Figure  2.6. 
The  experimental  data  are  plotted  as  the  experiment  progresses.  The  final  data  for  an 
experiment  conducted  on  a dummy  cell  are  shown  in  Figure  2.7  and  2.8.  Figure  2.7 
corresponds  to  an  experiment  conducted  with  settings  shown  in  Figure  2.5  (potentiostatic 
mode)  and  Figure  2.8  (galvanostatic  mode)  corresponds  to  the  settings  of  Figure  2.6.  A 
example  of  the  data  file  is  shown  in  Figure  2.9  and  2.10.  In  Figure  2.11  the  solid  line 
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represents  the  absolute  value  of  the  impedance  as  a function  of  frequency  for  the 
experiment  shown  in  Figure  2.5  and  the  dashed  line  represents  the  measuring  resistor 
used  to  collect  the  data.  The  measuring  resistor  switches  automatically  to  the  next  value 
when  the  change  in  the  absolute  value  of  the  impedance  shows  a significant  gradient.  The 
plot  of  the  absolute  value  of  impedance  and  the  measuring  resistor  as  a function  of 
frequency  for  the  data  shown  in  Figure  2.6  is  shown  in  Figure  2.12a.  The  measuring 
resistor  changes  twice  in  this  case.  The  first  change  in  the  measuring  resistor  occurs  in  the 
beginning  of  the  experiment  because  a suboptimal  measuring  resistor  was  chosen 
deliberately  at  the  start  of  the  experiment.  The  plot  of  the  perturbation  current  amplitude 
as  a function  of  frequency  is  shown  in  Figure  2.12b.  The  current  amplitude  changes  with 
the  change  in  the  absolute  value  of  impedance  to  maintain  linearity. 

2.4  Conclusions 

In  this  chapter  the  theory  of  impedance  measurements  was  summarized.  A brief 
description  of  the  instrumentation  and  the  principles  behind  their  working  was  outlined.  A 
menu  based  impedance  software  was  developed  in-house  to  conduct  impedance 
measurements.  The  unique  features  of  the  software  are  its  ability  to  choose  a measuring 
resistor  automatically  and  the  algorithm  to  make  galvanostatic  measurements. 
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Table  2.1: 

fhe  measuring  resistors  availa 

Die  on  SI1286|T81. 

Standard 

Resistor 

0.1 

1 

10 

100 

lk 

10k 

100k 

(Q) 

Full-scale 
Current 
(at  200  mv) 

2 

200m 

20m 

2m 

200g 

20  g 

±(A) 

Table  2.2:  The  measuring  resistors  available  on  PAR  273A[19]. 


Standard 

Resistor 

1 

10 

100 

lk 

10k 

100k 

1M 

10M 

(«) 

Full-scale 
Current 
(at  IV) 

1 

100m 

10m 

lm 

100(r 

10g 

lOOn 

±(A) 
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Figure  2. 1 : Small-signal  analysis  of  an  electrochemical  non-linear  system  [4]  (reproduced 
with  permission  from  the  author). 
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Figure  2.2:  Principle  of  correlation  analysis  [14]  (reproduced  with  permission  from 
author). 
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Figure  2.3:  a)  2-terminal  electrochemical  cell,  b)  3-terminal  electrochemical  cell,  c)  4- 
terminal  electrochemical  cell  [11]  (reproduced  with  permission  from  the  author). 
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Figure  2.4:  Schematic  of  the  experimental  setup. 
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Figure  2.6:  Control  front  panel  of  the  impedance  software  (Galvano static  Mode). 


Figure  2.7:  Data  collection  front  panel  for  Figure  2.5. 
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Figure  2.8:  Data  collection  front  panel  for  Figure  2.6. 
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-0.21561 

6500 

10.097 
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1000 
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1000 
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10 
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Figure  2.9:  Typical  data  file  for  an  experiment  with  settings  from  Figure  2.5. 
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0.092007 

109.4786 

-8.19745 

1.9953 

109.24 

-12.506 

0 

10000 

0.091342 

109.9535 

-6.5309 

1.5849 

109.78 

-10.124 

0 

10000 

0.090948 

110.2458 

-5.26896 

1.259 

110.13 

-7.9989 

0 

10000 

0.090706 

110.4201 

-4.15418 

1 

110.48 

-6.3804 

0 

10000 

0.090563 

110.6641 

-3.30525 

Figure  2. 10:  Typical  data  file  for  an  experiment  with  settings  from  Figure  2.6. 
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Figure  2.11:  The  solid  line  represents  the  absolute  value  of  the  impedance  for  the 
experiment  shown  in  Figure  2.5.  The  dashed  line  represents  the  measuring  resistor. 
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Figure  2.12:  a)  The  solid  line  represents  the  absolute  value  of  the  impedance  for  the 
experiment  shown  in  Figure  2.6.  The  dashed  line  represents  the  measuring  resistor,  b)  The 
circles  represent  the  perturbation  current  amplitude  as  a function  of  frequency. 


CHAPTER  3 


MEASUREMENT  MODELS:  DEMONSTRATION  OF 

APPLICABILITY 

3.1  Introduction 

Impedance  spectroscopy  has  become  a powerful  tool  for  investigating  the 
properties  of  electronic  materials  and  electrochemical  systems  [4,7],  The  response  of  a 
system  to  sinusoidal  perturbation  can  be  used  to  calculate  the  impedance  as  a function  of 
the  frequency  of  the  perturbation.  Interpretation  of  electrochemical  impedance  spectra 
requires  selection  of  an  appropriate  model  which  is  regressed  to  the  data.  The  most 
commonly  used  model  is  an  electrical  circuit  analogue  consisting  of  resistors,  capacitors, 
inductors,  and  specialized  distributed  elements.  Development  of  detailed  physico- 
electrochemical  models  is  becoming  more  common,  but  is  still  less  often  pursued. 
Selection  of  an  appropriate  model  can  be  time  consuming,  particularly  if  a detailed 
physico-electrochemical  model  is  desired.  Some  of  the  questions  that  must  be  addressed 
in  developing  an  appropriate  model  are: 

1.  Was  the  experimental  system  stationary?  One  distinction  between  electrochemical  and 
electrical  systems  is  that  electrochemical  processes  are  typically  dynamic,  i.e.  the  rate 
of  electrochemical  reactions  can  change  in  response  to  dissolution  of  the  electrode, 
deposition  of  salt  or  oxide  films,  and  adsorption  or  absorption  of  reacting  species.  If 
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the  system  changed  significantly  during  the  time  in  which  data  were  collected,  this 
nonstationary  behavior  should  be  incorporated  into  the  model. 

2.  How  many  physical  processes  can  be  resolved  from  the  data?  Electrochemical  systems 
rarely  display  well  resolved  peaks  in  the  impedance  spectra.  Since  the  ultimate  goal  of 
interpreting  impedance  spectra  is  to  gain  insight  into  the  physical  processes  that 
control  the  system,  the  model  developed  must  incorporate  time  constants  for  the 
relevant  physical  processes.  Another  way  to  phrase  this  question  is  to  ask  whether  the 
lack  of  complete  agreement  between  the  data  and  the  model  should  be  attributed  to 
experimental  error,  to  nonstationary  processes,  or  to  use  of  an  inadequate  or 
incomplete  model. 


3.2  Measurement  Models 


Models  can  be  classified  as  being  one  of  two  types.  Process  models  are  used  to 
predict  the  response  of  a system  from  physical  phenomena  that  are  hypothesized  to  be 
important.  Regression  of  process  models  to  data  allows  identification  of  physical 
parameters  based  upon  the  original  hypothesis.  In  contrast,  measurement  models  are  built 
by  sequential  regression  of  line  shapes  to  the  data.  The  model  can  be  used  to  identify 
characteristics  of  the  data  set  that  could  facilitate  selection  of  an  appropriate  process 
model. 


The  measurement  model  selected  for  this  work  was  the  Voigt  model 


Z(o))  = Z0  + 


y_A__ 

* (1  + jTka>) 


(3.1) 
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An  electrical  circuit  corresponding  to  the  Voigt  model  is  presented  in  Figure  3.1.  Use  of 
the  Voigt  model  was  motivated  by  the  routine  use  of  measurement  models  in  the 
deconvolution  of  optical  spectra  [22],  To  be  useful,  the  measurement  model  must  contain 
line  shapes  that  correspond  to  physical  processes.  Since  stationary  optical  (and  impedance) 
spectra  satisfy  the  Kramers-Kronig  relations,  the  model  must  also  be  consistent  with  the 
Kramers-Kronig  relations.  The  impedance  response  of  a single  Voigt  element  is  usually 
attributed  to  the  electrical  response  of  a linearized  electrochemical  reaction  [7].  The  term 
Z0  corresponds  to  the  ohmic  resistance,  the  time  constant  it  for  element  k is  equivalent  to 
RkCic  in  the  Voigt  model,  and  At  is  equivalent  to  Rjc.  Equation  (3.1)  is  consistent  with  the 
Kramers-Kronig  relations. 

The  tenets  underlying  the  use  of  measurement  models  for  impedance  spectroscopy 
are: 

1.  By  including  a sufficient  number  of  terms,  a general  measurement  model  based  on 
equation  (3.1)  can  fit  impedance  data  for  typical  stationary  electrochemical  systems. 
The  object  of  this  chapter  is  to  demonstrate  that  this  postulate  holds  for 
electrochemical  systems  that  exhibit  the  influence  of  mass  transfer  (through  Warburg 
elements),  pseudo-capacitive  behavior,  frequency  dispersion  (through  constant  phase 
elements),  and  kinetic  control. 

2.  Because  the  measurement  model  does  fit  stationary  impedance  data,  an  inability  to 
predict  the  imaginary  part  of  the  impedance  spectrum  with  parameters  obtained  by 
fitting  the  real  part  of  the  spectrum  (or,  conversely,  to  predict  the  real  part  of  the 
impedance  spectrum  with  parameters  obtained  by  fitting  the  imaginary  part  of  the 
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spectrum)  can  be  attributed  to  failure  of  the  data  to  conform  to  the  assumptions  of  the 
Kramers-Kronig  relations  rather  than  to  failure  of  the  model.  Thus,  as  discussed  in 
references  [23]  and  [24]  the  measurement  model  could  be  used  to  assess  the 
consistency  of  impedance  data  with  the  Kramers-Kronig  relations  without  integration. 
3.  Another  result  of  the  adequacy  of  the  measurement  model  is  that  the  model  can  be 
used  to  identify  the  frequency-dependent  error  structure  of  impedance  spectra.  The 
error  structure  can  then  be  used  to  weight  the  data  during  regression  and  to  provide  a 
means  of  deciding  whether  a given  regression  provided  a “good  fit”  [23,24], 

The  use  of  measurement  models  is  superior  to  the  use  of  polynomial  fitting  because  fewer 
parameters  are  needed  to  model  complex  behavior,  and  because  the  measurement  model 
satisfies  the  Kramers-Kronig  relations  implicitly.  Experimental  data  can,  therefore,  be 
checked  for  consistency  with  the  Kramers-Kronig  relations  without  actually  integrating  the 
equations  over  frequency.  The  use  of  measurement  models  does  not  require  extrapolation 
of  the  experimental  data  set;  therefore,  inaccuracies  associated  with  an  incomplete 
frequency  spectrum  are  resolved.  For  the  application  to  a preliminary  screening  of  the 
data,  the  use  of  measurement  models  is  superior  to  the  use  of  more  specific  electrical 
circuit  analogues  because  one  can  determine  whether  the  residual  errors  are  due  to  an 
inadequate  model,  to  failure  of  data  to  conform  to  the  Kramers-Kronig  assumptions,  or  to 
experimental  noise. 

3.3  Approach  for  Demonstration  of  Applicability 
The  usefulness  of  the  measurement  model  depends  on  the  extent  to  which  the  line 
shape  for  each  element  corresponds  to  that  generated  by  physical  phenomena.  For 
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example,  the  utility  of  measurement  models  in  the  field  of  optics  is  due  to  the  accuracy  of 
the  Lorentzian,  Drude,  and  Debye  models  of  optical  excitation  [26-28],  The  claims  made 
here  for  the  measurement  models  require  that  the  line  shape  used  be  adequate  to  represent 
typical  impedance  data.  While  this  approach  has  been  applied  directly  to  experimental 
data,  [24,25]  the  intent  here  is  to  demonstrate  applicability  to  published  data  for  some 
rather  complex  systems.  The  authors  cited  have  used  electrical  circuit  analogues  to  model 
their  results,  and,  in  this  work,  the  measurement  models  were  applied  to  synthetic  data 
obtained  from  the  authors'  models.  The  models  were  selected  for  this  work  because  they 
encompass  phenomena  typically  seen  in  electrochemical  impedance  spectra  (i.e.,  diffusion, 
pseudo-capacitive  behavior,  and  frequency  dispersion).  The  circuit  analogues  are 
presented  in  Figure  3.2  and  the  associated  parameter  values  are  presented  in  Table  3.1. 

The  “synthetic  data”  to  which  the  measurement  model  was  applied  were  calculated 
in  double  precision  from  the  models  presented  in  Figure  3.2  and  Table  3.1.  The  data  were, 
therefore,  inherently  consistent  with  the  Kramers-Kronig  relations,  and  the  only  “noise”  in 
the  data  was  caused  by  round-off  error  in  the  calculation.  Use  of  synthetic  data  provides  a 
good  test  of  the  measurement  model  because  the  errors  associated  with  experimental  noise 
and  nonstationary  behavior  are  avoided.  A poor  fit  of  the  measurement  model  to  synthetic 
data  can,  therefore,  be  attributed  only  to  inadequacies  of  the  model.  Two  nonlinear 
regression  packages,  CNLS  by  J.  R.  Macdonald  [8]  and  a proprietary  program  written  by 
L.  Garcia-Rubio  [34],  were  used  to  fit  the  measurement  model  to  the  synthetic  data,  and 
the  fits  obtained  were  identical.  The  data  were  normalized  such  that  each  point  had  equal 
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weight  in  the  regression.  This  normalization  was  used  to  avoid  excessive  weighting  on  the 
low  frequency  data  which  typically  have  much  larger  values  for  the  impedance. 

3.4  Results  and  Discussion 

The  sequential  manner  in  which  a measurement  model  is  constructed  is  demonstrated  in 
the  first  subsection.  The  application  of  the  measurement  model  is  further  demonstrated  for 
systems  influenced  by  diffiisional  processes  and  inductive  behavior.  The  applicability  of 
the  measurement  model  is  also  demonstrated  for  a system  that  was  modeled  by  a constant 
phase  element  and  for  a high  impedance  system.  The  interpretation  of  parameter  values  is 
discussed  in  the  final  subsection. 

3 .4.1  Sequential  Construction  of  a Model:  Circuit  1 

In  accordance  with  equation  (3.1),  a measurement  model  is  constructed  by 
sequentially  adding  k Voigt  elements  with  parameters  Ak  and  xk  until  the  fit  is  no  longer 
improved  by  addition  of  yet  another  element.  This  procedure  is  illustrated  for  synthetic 
data  obtained  from  the  electrical  circuit  model  proposed  for  evolution  of  hydrogen  on  a 
LaNis  ingot  electrode  [24],  A diffiisional  resistance  was  observed  in  this  system  as  a result 
of  the  incorporation  of  hydrogen  into  the  metal.  The  electrical  circuit  model  presented  as 
Circuit  1 in  Figure  3.2  accounted  for  mass  transfer  of  hydrogen  with  a Warburg  element, 
and  three  time  constants  are  evident  in  the  circuit  analogue.  The  complex  Warburg 
impedance,  given  by 

-7  _ D tanh(7/c^©) 

(v;v®) 


(3.2) 
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can  be  derived  from  the  solution  of  V2c(  = 0 for  a film  of  thickness  6 [4,7],  The  time 

constant  can  be  expressed  in  terms  of  diffusivity  and  film  thickness  as  xw  = . The 

Warburg  element  is  sometimes  called  a distributed  element  because  it  can  be  modeled  by  a 
distribution  of  relaxation  times  [7], 

The  manner  in  which  the  measurement  model  behaves  upon  successive 
incorporation  of  line  shapes  is  demonstrated  in  Figures  3.3  through  3.7.  The  difference 
between  the  simulated  data  and  a measurement  model  incorporating  one  Voigt  element  is 
presented  in  Figure  3.3.  The  residual  sum  of  squares  obtained  by  the  fit  was  37.15.  The 
residual  sum  of  squares  was  improved  to  0.5588  by  addition  of  a second  line  shape  as  seen 
in  Figure  3.4.  Inclusion  of  a third  line  shape,  shown  in  Figure  3.5,  resulted  in  a three  order 
of  magnitude  decrease  in  the  residual  sum  of  squares.  The  improvement  of  the  fit 
achieved  by  addition  of  a fourth  element  (Figure  3.6)  is  not  readily  seen  in  a comparison  of 
Figures  3.5  and  3.6.  The  inclusion  of  a fourth  element  was,  however,  considered  to  be 
statistically  significant  because  the  residual  sum  of  squares  was  improved  by  an  additional 
three  orders  of  magnitude.  A significant  improvement  in  the  fit  was  also  observed  for  five 
Voigt  elements  (Figure  3.7),  and  the  resultant  residual  sum  of  squares  was  4.915  xlO'10 

The  residual  sum  of  squares  obtained  by  fitting  the  measurement  model  to  the 
synthetic  data  is  presented  as  a function  of  the  number  of  Voigt  elements  in  Figure  3.8. 
Figure  3.8  can  be  used  to  show  that  a maximum  of  five  Voigt  elements  can  be  resolved 
from  the  synthetic  data  because  the  sixth  Voigt  element  did  not  improve  the  quality  of  the 
fit.  Figure  3.8  also  shows  that,  while  visual  inspection  of  Figure  3.5  suggests  that  the  fit 
obtained  with  three  elements  was  acceptable  (resulting  in  residual  errors  less  than  0.1 
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percent),  addition  of  the  fourth  and  fifth  element  provided  significant  improvement  in  the 
residual  error.  The  residual  error  obtained  by  regressing  the  five-time-constant  Voigt 
model  to  the  data  is  presented  as  a function  of  frequency  in  Figure  3.9.  The  average 
residual  error  was  less  than  2 x 104  percent.  Several  points  must  be  made  here.  The 
measurement  model  based  on  a summation  of  Voigt  elements  provided  an  excellent  fit,  in 
spite  of  the  fact  that  the  circuit  included  a Warburg  element  which  provides  an  impedance 
response  that  differs  significantly  from  that  obtained  for  a resistor  and  capacitor  in  parallel. 
The  measurement  model  fit  the  synthetic  data  to  within  7 significant  figures;  therefore,  the 
measurement  model  can  be  considered  to  be  statistically  adequate  to  fit  the  corresponding 
experimental  impedance  data  for  which  only  3 or  4 significant  figures  can  be  expected. 

To  explore  the  potential  influence  of  experimental  error,  random  noise  with  a 
maximum  amplitude  of  five  percent  was  added  to  the  results  of  Circuit  1 (and  reported  as 
Circuit  lb).  The  optimal  regression  of  the  measurement  model  yielded  only  three  Voigt 
circuit  elements  as  shown  in  Figure  10.  As  shown  in  Figure  3.8,  only  three  elements  could 
be  resolved  because  the  contribution  of  the  fourth  circuit  element  was  insignificant  as 
compared  to  the  noise  of  the  data.  The  normalized  sum  of  squares  obtained  by  regression 
of  the  measurement  model  was  slightly  lower  than  the  noise  level  which  is  shown  as  a 
dashed  line  in  Figure  3.9.  The  residual  errors  are  presented  as  a function  of  frequency  in 
Figure  3.11  and  are  randomly  distributed  with  respect  to  frequency.  The  parameter  values 
obtained  through  regression  of  the  measurement  model  to  the  LaNi5  system  are  presented 
in  Table  3.2.  The  parameter  values  obtained  for  the  first  three  elements  changed  slightly 
(about  5 percent)  to  compensate  for  the  noise  and  the  lack  of  the  fourth  and  fifth  elements. 
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As  expected,  only  a rough  equivalency  is  seen  between  the  fit  parameters  and  the 
parameters  of  the  original  circuit.  For  example,  the  largest  time  constant  obtained  in  the 
fit  (Table  3.2)  was  0.387s;  whereas,  the  largest  time  constant  in  the  original  circuit  (Table 
3.1)  had  a value  of  0.7055s.  Good  agreement  was  seen  between  the  Ohmic  resistance  of 
the  original  circuit  R,  and  the  high  frequency  limit  of  the  measurement  model  Z0. 

3.4.2  Application  to  Diffusion  Systems:  Circuit  2 

In  general,  the  number  of  Voigt  elements  needed  to  model  a system  which  is 
affected  by  diffusion  exceeds  the  number  of  time  constants  implicit  in  a model  that 
incorporates  a Warburg  element.  While  the  contribution  of  the  additional  Voigt  elements 
can,  as  shown  above,  be  insignificant  in  comparison  to  experimental  noise,  the  discrepancy 
remains  because  a single  Voigt  element  simply  cannot  account  for  the  frequency 
dependence  associated  with  mass  transfer.  If  a measurement  model  is  to  be  used  to 
identify  the  correct  number  of  time  constants  for  a given  system,  additional  model 
elements  must  be  incorporated  in  much  the  same  way  as  multiple  models  (Lorentzian, 
Drude,  or  Debye)  can  be  used  to  develop  measurement  models  for  optical  dispersion.  The 
number  of  model  elements  used  is,  of  course,  irrelevant  if  the  measurement  model  is  to  be 
used  simply  to  check  the  consistency  of  data  to  the  Kramers-Kronig  assumptions. 

The  application  of  the  Voigt  measurement  model  is  presented  in  Figure  3.12  for 
the  corrosion  of  carbon  steel  in  3wt%  NaCl  solution  near  the  corrosion  potential  [29],  The 
corresponding  electrical  circuit  model  (presented  as  Circuit  2 in  Figure  3.2)  included  a 
Warburg  diffusion  element  to  account  for  the  mass  transfer  resistance  of  the  solution. 
Eight  Voigt  elements  were  needed,  even  though  only  three  time  constants  are  present  in 
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circuit  2.  The  results  of  this  regression  do  show  that  the  measurement  model  can  provide 
an  adequate  representation  of  experimental  impedance  data  that  can  be  modeled  in  terms 
of  a Warburg  impedance.  Similar  agreement  was  obtained  for  direct  analysis  of 
experimental  data  for  systems  influenced  by  mass  [24,25], 

3.4.3  Application  to  Inductive  Systems:  Circuit  3 

The  appearance  of  a positive  imaginary  component  of  the  impedance  has  been 
attributed  at  high  frequencies  to  the  stray  capacity  of  the  current  measuring  resistor  and  at 
low  frequencies  to  adsorption  phenomena  at  the  electrode  surface  [7],  Transport  based 
models  for  this  behavior  are  rare,  and  most  electrical  circuit  models  account  for  this 
behavior  by  incorporating  either  an  inductor  or  a capacitor  with  a negative  value  for  the 
capacitance.  The  measurement  model  was  applied  to  pseudo-capacitive  data  reported  by 
Lorenz  and  Mansfeld  [30]  for  corrosion  of  iron  in  0.5M  H2SO4.  Their  model,  shown  as 
Circuit  3 in  Figure  3.2,  incorporated  an  inductor  to  fit  the  pseudo-capacitive  response  at 
low  frequencies.  The  Voigt  model  can  account  for  inductive  behavior  if  the  magnitude  Ak 
is  allowed  to  have  negative  sign.  The  regression  to  Circuit  3 is  given  in  Figure  3.13.  The 
measurement  model  indicates  that  three  time  constants  can  be  resolved  from  the  synthetic 
data,  and  this  result  is  fully  consistent  with  the  circuit.  Complete  agreement  was  obtained 
because  the  circuit  is  composed  of  only  passive  elements  (one  inductor,  three  resistors, 
and  two  capacitors).  The  residual  sum  of  squares  for  this  fit  was  found  to  be  3.221  x 10'14. 
The  measurement  model  provides  a good  fit  to  data  exhibiting  an  inductive  response. 
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3.4.4  Application  to  Constant  Phase  Elements:  Circuit  4 


The  broadening  or  depression  of  a semicircle  in  the  impedance  plane  plot  can  be 
attributed  to  a variety  of  physical  phenomena.  For  example,  models  have  been  developed 
that  associate  depressed  semicircles  to  multiple  or  coupled  reaction  sequences,  to 
roughening  of  the  electrode,  and  to  frequency-dependent  ohmic  resistances  caused  by  a 
non-uniform  charging  of  the  double  layer  [7],  This  type  of  impedance  response  cannot  be 
simulated  by  simply  incorporating  a Warburg  term  in  the  model;  therefore,  generalized 
impedance  elements  termed  “constant-phase  elements  (CPE)”  are  often  incorporated  into 
the  model.  The  measurement  model  was  applied  to  the  results  of  a circuit  proposed  by 
Kendig  et  al.  [31]  for  a model  pit  electrode.  The  electrical  circuit  is  presented  as  Circuit  4 
in  Figure  3.2.  The  mathematical  formulation  used  for  the  constant  phase  element  was 


R „ 


'CPE 


(1  + (7'xCP£co)p 


(3.3) 


whereP  has  value  between  0 and  1.  As  seen  in  Figure  3.14,  with  addition  of  a sufficient 
number  of  elements,  the  Voigt  measurement  model  can  provide  a good  representation  of  a 
constant  phase  element.  Eight  Voigt  elements  were  used  in  the  regression;  whereas,  only 
two  time  constants  are  evident  in  the  circuit  itself.  As  shown  earlier,  experimental  noise 
will  reduce  the  number  of  Voigt  elements  that  can  be  resolved. 

3.4.5  Application  to  Coated  Metals:  Circuit  5 

The  measurement  model  presented  here  is  most  effective  in  modeling  systems  that 
can  be  modeled  otherwise  with  passive  elements  such  as  resistors,  capacitors,  and 
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inductors.  It  is  perhaps  no  surprise  that,  with  only  two  elements,  the  measurement  model 
can  provide  a good  representation  of  a circuit  developed  for  painted  metals  that  does  not 
include  diffusional  or  constant  phase  elements.  The  electrical  circuit  is  presented  as  Circuit 
5 in  Figure  3.2.  The  fit  to  synthetic  data  obtained  from  a model  for  painted  steel  is 
presented  in  Figure  3.15.  The  normalized  sum  of  squares  for  this  fit  was  found  to  be 
2.055  xlO'14. 

3.4.6  Interpretation  of  Parameter  Values 

It  is  not  intended  that  use  of  measurement  models  replace  the  development  of 
physico-electrochemical  models  or  circuit  analogues  specific  to  a given  system.  Instead, 
the  use  of  a measurement  model  is  intended  to  guide  selection  and  development  of  more 
specific  models.  The  only  generalization  that  can  be  made  for  electrochemical  systems  is 
that  the  high  frequency  asymptote  Z0  corresponds  to  the  ohmic  resistance  R,.  It  may  be 
possible  to  derive  a relationship  between  the  parameters  Tk  and  Ak  and  physical 
parameters,  but  such  a relationship  must  be  developed  for  each  system. 

3.5  Conclusions 

A measurement  model  based  on  the  Voigt  circuit  was  shown  in  this  work  to 
provide  a statistically  adequate  fit  to  models  that  represent  typical  impedance  data.  The 
measurement  model  can  be  used  to  check  the  consistency  of  electrochemical  impedance 
data  to  the  Kramers-Kronig  relations  because  the  model  itself  satisfies  the  Kramers-Kronig 
relations  and  because  the  model  provides  a statistically  adequate  fit  to  electrochemical 
impedance  data.  The  approach  taken  would  be  to  fit  the  model  to  either  the  real  or  the 
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imaginary  component  of  the  spectrum  and  to  predict  the  other  component  from  the 
regressed  model  parameters.  The  applicability  of  the  model  established  here  allows  use  of 
the  standard  deviation  of  model  parameters  to  establish  a confidence  window  for  the 
predicted  component.  Thus,  a sound  statistical  basis  can  be  established  for  rejection  of 
data  as  being  inconsistent  with  a stationary  process.  This  approach  will  be  presented  in  a 
chapter  5. 

The  measurement  model  provided  a one-to-one  identification  of  the  time  constants 
for  a class  of  processes  that  are  described  in  terms  of  passive  elements.  Thus,  the 
measurement  model  may  allow  correct  identification  of  the  number  of  physical  processes 
that  can  be  discerned  from  the  impedance  response  of  systems  that  are  unaffected  by  mass 
transfer  or  by  the  frequency  dispersion  associated  with  nonuniform  current/potential 
distributions.  Solid  state  systems  may,  for  example,  be  suited  for  this  type  of  analysis.  The 
measurement  model  presented  here  cannot,  however,  be  generally  used  for  identification 
of  the  largest  number  of  time  constants  discernible  in  a given  spectrum  because  several 
Voigt  elements  were  needed  to  model  a distributed  element. 
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Table  3.1:  Val 

ues  for  circuit  parameters 

Circuit: 

l(aand  b)  (16) 

2(28) 

3(29) 

4(30) 

5(31) 

23.22 

- 

0.5 

1000 

Rj.fi 

418.3 

150 

100 

15 

2 x 105 

C,,uF 

18.82 

1000 

25 

128 

1 x 103 

Rz.O 

695.2 

150 

200 

- 

4x  105 

C2,mF 

92.58 

- 

50 

- 

1 x 10'2 

r3,o 

- 

- 

100 

- 

- 

L,H 

- 

- 

100 

- 

_ 

Rw,  fi 

848.4 

500 

- 

- 

. 

Tw,  S 

0.7055 

5 

- 

- 

. 

Rp,  fi 

- 

- 

- 

1135 

"Tcpe,  s 

- 

- 

- 

10 

p 

- 

- 

- 

0.79 
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Table  3.2:  The  results  of  the  Regression  of  the  Voigt  Measurement  model  (Fig.  3.1)  to 
the  circuits  shown  in  Fig  3.2. 


Circuit 

la 

(Warburg) 

lb 

(Warburg) 

2 

(Warburg) 

3 

Inductor 

4 

CPE 

5 

Zo,Q 

23.22 

23.77 

0 

0 

0.50036 

1000 

A,,Q 

266.48 

289.7 

10.60 

58.695 

274.25 

436900 

Tj,  S 

0.006404 

0.006569 

0.7954 

0.002932 

3.339 

0.00442 

a2,q 

1163 

1201 

28.04 

-102.05 

13.1 

163100 

t2,  s 

0.388 

0.3870 

4.101 

343.35 

0.001823 

0.000181 

A3,  Q 

494.0 

528.8 

22.83 

0.01709 

3.323 

t3,  s 

0.0760 

0.0731 

0.1131 

0.01164 

A,), 

25.61 

51.50 

11.29 

t4,  s 

0.02723 

0.07162 

0.09767 

As,  £2 

9.126 

8.004 

49.91 

ts,  s 

0.00854 

0.2877 

0.6158 

A6,Q 

0.8880 

162.62 

Tg,  S 

0.03041 

57.67 

A7,n 

0.01248 

605.75 

x7,  s 

0.006603 

12.15 

As,  a 

0.00003352 

26.01 

tg,  s 

0.0004966 

492.9 

1.  Circuit  1:  hydrogen  evolution  at  a LaNi5  electrode  in  30wt%  KOH  solution  at  an 
applied  potential  of -1.  IV  (Hg/HgO)  [24]  Circuit  lb  has  5%  random  noise  added. 

2.  Circuit  2:  corrosion  of  carbon  steel  in  3wt%  NaCl  solution  near  the  corrosion  potential 
[29J. 

3 . Circuit  3 : corrosion  of  iron  in  0. 5M  H2SO4  at  the  corrosion  potential  [30], 

4.  Circuit  4:  corrosion  of  a model  pit  electrode  in  0.5M  NaCl  [31], 

5.  Circuit  5:  corrosion  of  a painted  metal  [32], 
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Figure  3.1:  The  Voigt  electrical  circuit  analogue  corresponding  to  the  measurement  model 
used  in  this  work  (equation  (3.1)). 
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Circuit  1: 


Circuit  2: 


Circuit  3: 


h 


Circuit  4: 


Circuit  5: 


Figure  3.2:  Electrical  circuit  analogues  for:  1)  hydrogen  evolution  on  LaNi5  [24];  2) 
corrosion  of  carbon  steel  in  3wt%  NaCl  solution  near  the  corrosion  potential  [29];  3) 
corrosion  of  iron  in  0.5M  H2SO4  at  the  corrosion  potential  [30];  4)  corrosion  of  a model 
pit  electrode  in  0.5M  NaCl  [3 1];  and  5)  corrosion  of  a painted  metal  [32], 
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Figure  3.3:  Regression  of  a measurement  model  with  one  Voigt  element  to  synthetic  data 
(o)  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNi5  [24]).  Solid  lines  in  all  figures 
are  the  results  of  the  model  fit.  The  residual  sum  of  squares  for  the  fit  was  37. 15. 
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Frequency,  Hj 


Figure  3.4:  Regression  of  a measurement  model  with  two  Voigt  elements  to  synthetic  data 
(o)  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNi5  [24]).  The  residual  sum  of 
squares  for  the  fit  was  0.5588.  Solid  lines  in  all  figures  are  the  results  of  the  model  fit. 
The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its  deconvolution  into 
the  contribution  of  each  Voigt  element. 
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Figure  3.5:  Regression  of  a measurement  model  with  three  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNi5  [24]).  The  residual  sum  of 
squares  for  the  fit  was  2.112  x 1CT4  Solid  lines  in  all  figures  are  the  results  of  the  model 
fit.  The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its  deconvolution 
into  the  contribution  of  each  Voigt  element. 
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Figure  3.6:  Regression  of  a measurement  model  with  four  Voigt  elements  to  synthetic  data 
(o)  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNij  [24]).  The  residual  sum  of 
squares  for  the  fit  was  1.772  xlO'7.  Solid  lines  in  all  figures  are  the  results  of  the  model  fit. 
The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its  deconvolution  into 
the  contribution  of  each  Voigt  element. 
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frequency,  Hz 

Figure  3.7:  Regression  of  a measurement  model  with  five  Voigt  elements  to  synthetic  data 
(o)  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNi5  [24]).  The  residual  sum  of 
squares  for  the  fit  was  4.915  x 10'10.  Solid  lines  in  all  figures  are  the  results  of  the  model 
fit.  The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its  deconvolution 
into  the  contribution  of  each  Voigt  element. 


Sum  of  Squores 
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Figure  3.8:  Normalized  sum  of  squares  for  the  regression  of  a measurement  model  to 
synthetic  data  obtained  from  Circuit  1 (hydrogen  evolution  on  LaNi5  [24])  as  a function  of 
the  number  of  Voigt  time  constants  employed  in  the  model:  (o)  regression  to  synthetic 
data;  (A)  regression  to  synthetic  data  with  random  noise  added;  and  (dashed  line) 
normalized  noise  level. 
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Figure  3.9:  Residual  errors  (o)  as  a function  of  frequency  for  the  regression  of  a 
measurement  model  with  five  Voigt  elements  to  synthetic  data  obtained  from  Circuit  1 
(hydrogen  evolution  on  LaNis  [24]). 
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Figure  3.10:  Regression  of  a measurement  model  with  three  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  2 with  added  random  noise(5%  amplitude).  The  residual 
sum  of  squares  for  the  fit  was  2. 424x1  O'2.  Solid  lines  in  all  figures  are  the  results  of  the 
model  fit.  The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its 
deconvolution  into  the  contribution  of  each  Voigt  element. 
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Figure  3.11:  Residual  errors  as  a function  of  frequency  for  the  regression  of  a 
measurement  model  with  three  Voigt  elements  to  synthetic  data  (o)  obtained  from  Circuit 
2 with  added  random  noise  (5%  maximum  amplitude). 
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Figure  3.12:  Regression  of  a measurement  model  with  four  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  2 (corrosion  of  carbon  steel  in  3 wt%  NaCl  solution  near 
the  corrosion  potential  [29]).  The  residual  sum  of  squares  for  the  fit  was  3.763X1CT6.  Solid 
lines  in  all  figures  are  the  results  of  the  model  fit.  The  lines  presented  in  the  lower  two 
figures  are  the  model  fit  and  its  deconvolution  into  the  contribution  of  each  Voigt  element. 
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Figure  3.13:  Regression  of  a measurement  model  with  three  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  3 (corrosion  of  iron  in  0.5M  H2SO4  at  the  corrosion 
potential  [30]).  The  residual  sum  of  squares  for  the  fit  was  3.221xl0'14.  Solid  lines  in  all 
figures  are  the  results  of  the  model  fit.  The  lines  presented  in  the  lower  two  figures  are  the 
model  fit  and  its  deconvolution  into  the  contribution  of  each  Voigt  element. 
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Figure  3.14:  Regression  of  a measurement  model  with  eight  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  4 (corrosion  of  a model  pit  electrode  in  0.5M  NaCl  [30]). 
The  residual  sum  of  squares  for  the  fit  was  1.49 lx  10'2.  Solid  lines  in  all  figures  are  the 
results  of  the  model  fit.  The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and 
its  deconvolution  into  the  contribution  of  each  Voigt  element. 
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Figure  3.15:  Regression  of  a measurement  model  with  two  Voigt  elements  to  synthetic 
data  (o)  obtained  from  Circuit  5 (corrosion  of  a painted  metal  [32]).  The  residual  sum  of 
squares  for  the  fit  was  2.055x  10'u.  Solid  lines  in  all  figures  are  the  results  of  the  model  fit. 
The  lines  presented  in  the  lower  two  figures  are  the  model  fit  and  its  deconvolution 


CHAPTER  4 


MEASUREMENT  MODELS:  DETERMINATION  OF  ERROR 

STRUCTURE 


In  the  previous  chapter,  it  was  shown  that  a measurement  model  based  on  Voigt 
circuit  elements  can  fit  typical  electrochemical  impedance  spectroscopy  data.  In  this 
chapter,  the  measurement  model  is  used  to  identify  the  frequency  dependent  error 
structure  of  impedance  data.  An  empirical  model  for  the  stochastic  component  of  the  error 
is  proposed.  It  is  shown  how  the  model  for  the  error  structure  can  be  used  to  guide  the 
design  of  experiments  and  used  to  weight  the  data  during  regression. 


4 . 1 Classification  of  Measurement  Errors 

The  residual  errors  (zresiduai)  that  arise  due  to  the  regression  of  a model  (Z)  to 
experimental  data  (Zexp)  can  be  of  two  types,  systematic  errors  (e systematic)  and  random  or 
stochastic  errors  (zst0chastic)  [34],  more  specifically. 


^exp  ^ + £ residual 


^ residual  ^ systematic  + ^ stochastic 


(4.1) 


The  systematic  errors  can  arise  due  to  the  lack  of  fit  of  the  model  (£/„/)  to  the  data  or  some 
kind  of  bias  (e bias)  in  the  experiment,  i.e., 


^ sytematic  ^ lof  ~^~^bias 


(4-2) 


71 


72 


In  electrochemical  impedance  spectroscopy  the  systematic  errors  ebias  can  arise  if 
the  system  is  non-stationary  (e„),  i.e.,  is  changing  during  the  course  of  an  experiment,  or 
due  to  instrumental  artifacts  (e,„s): 

e6ias=em+81>u  (4.3) 

It  is  generally  accepted  that  ens  are  observed  at  low  frequencies,  where  relatively  long 
periods  of  time  are  required  to  collect  data.  In  addition,  e ins  are  typically  considered  high- 
frequency  effects.  Most  electrochemical  systems  are  inherently  non-stationary  due  to  the 
change  in  the  electrode  surface  during  the  course  of  an  experiment  (the  surface  could 
change  due  to  corrosion,  electrodeposition,  film  formation  etc.).  In  contrast,  solid-state 
systems  can  be  assumed  to  be  stationary  as  a first  approximation.  Impedance  data  can  be 
corrupted  by  instrumental  artifacts  in  electrochemical  systems  and  solid  state  systems 
alike.  In  impedance  spectroscopy  the  presence  of  ebias  results  in  the  data  being 
inconsistent  with  the  Kramers-Kronig  relations.  The  issues  that  arise  out  of  these 
inconsistencies  will  be  discussed  in  the  next  chapter. 

The  presence  of  stochastic  errors  z.tochaaic  in  any  experimental  data  is  inevitable.  In 
this  chapter,  8 stochastic  are  discussed.  A technique  to  filter  out  8 stochastic  from  er«/dw  is 
outlined  and  an  empirical  model  for  est0chastic  is  proposed.  Impedance  data  can  range  over 
several  orders  of  magnitude,  which  results  in  the  data  being  strongly  heteroscedastic. 
Hence  a priori  knowledge  of  the  error  structure  of  impedance  data  leads  to  improvements 
in  the  regression  of  a model  to  the  data  by  using  proper  weighting. 
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4.2  Determination  of  Stochastic  Errors 


In  this  chapter,  unless  explicitly  indicated,  Z will  refer  to  a measurement  model 


based  on  Voigt  circuits.  Combination  of  equations  4. 1-4.3  results  in 


stochastic 


(4.4) 


In  the  preceding  chapter,  it  was  shown  that  a measurement  model  based  on  Voigt 
circuit  elements  can  adequately  represent  electrochemical  impedance  data.  As  a result,  by 
adding  a sufficient  number  of  Voigt  line  shapes,  and  by  using  optimal  weighting  (see 
section  4.4),  it  can  be  assumed  that  S/0/is  negligible  compared  to  the  stochastic  component 


accuracy  of  this  assumption. 

4.2.1  Stationary  Systems 

The  procedure  for  calculation  of  the  stochastic  error,  estochastic,  in  stationary  systems 
is  illustrated  by  an  example  of  impedance  data  for  n-GaAs  single  crystal  wafer  with  a Ti 
Schottky  contact  and  a Au/Ge/Ni  ohmic  contact  [35] . Experimental  data  for  five  replicate 
experiments  conducted  at  320K  are  shown  in  Figure  4.1.  For  stationary  systems  (see 
equation  4.4),  era  is  equal  to  zero.  Therefore  equation  4.4  reduces  to, 


replicate  data.  Under  the  assumption  that  the  stochastic  errors  follow  a normal  distribution 
with  zero  mean  and  variance  o2  [36], 


of  the  noise.  The  accuracy  of  the  algorithms  outlined  in  this  section  will  depend  on  the 


(4.5) 


A 

where  Z is  calculated  by  regression  of  a measurement  model  to  the  combined  set  of 


^ stochastic  ~ O') 


(4.6) 
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From  equation  4.6  the  mean  of  the  stochastic  error  is  equal  to  f , and,  since  for  replicate 
experiments  the  errors  due  to  instrumental  artifact  can  be  assumed  to  be  constant  from 
one  experiment  to  another,  the  frequency-dependent  standard  deviation,  a can  be 
calculated  by  taking  the  standard  deviation  of  the  residual  error  erejWufl/  at  each  frequency, 


^.2  _ y.  residual,. k ^ residual ) 

~ N 


(4.7) 


where  ct2  is  the  calculated  variance  of  the  residual  errors  and  N is  the  number  of  data 
points. 

The  results  of  the  calculation  of  ct  for  the  GaAs  data  in  Figure  4. 1 are  shown  in 
Figure  4.2.  The  circles  in  Figure  4.2  represent  the  real  part  of  the  standard  deviation,  ct 
and  the  triangles  represent  the  imaginary  part.  Examination  of  the  figure  shows  that  the 
real  and  imaginary  part  of  ct  are  equal.  This  is  an  important  result,  and  its  significance  will 
be  discussed  in  the  later  sections. 


4.2.2  Non-Stationary  Systems 


For  non-stationary  systems  the  error  ens  is  not  equal  to  zero.  In  fact,  since  the 
system  evolves  with  time,  e„  changes  from  one  experiment  to  another  as  a set  of  replicate 
(or  consecutive)  experiments  are  performed.  Hence,  equation  4.7  leads  to  inaccurate 
estimates  for  ct  because  the  contribution  of  to  the  standard  deviation  cannot  be  ignored. 

One  can  identify  three  time  scales  over  which  an  electrochemical  system  can 
change  during  the  course  of  an  experiment.  Usually  several  consecutive  impedance 
experiments  are  conducted  where  data  is  collected  at  several  frequencies  for  each  scan. 
The  smallest  time  scale  is  the  time  to  collect  a data  point  at  one  frequency.  The  time 
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required  to  collect  one  complete  impedance  scan  constitutes  the  next  time  scale,  and  the 
third  time  scale  is  the  total  time  elapsed  from  the  start  of  the  first  experiment  to  the  end  of 
the  last  experiment. 

If  the  system  is  evolving  very  rapidly,  then  the  change  can  occur  during  the  time  in 
which  one  data  point  is  collected.  The  likelihood  of  the  system  changing  during  the 
collection  of  single  data  point  is  high  at  low  frequencies  where  the  time  required  to  make 
one  measurement  is  high  (inverse  of  the  frequency).  For  slower  changing  systems,  the 
change  might  not  be  apparent  at  each  frequency,  but  significant  change  can  occur  between 
the  start  and  end  of  one  complete  frequency  scan.  These  kinds  of  non-stationarities  result 
in  the  data  being  inconsistent  with  the  Kramers-Kronig  relations.  The  issues  arising  out  of 
these  inconsistencies  will  be  discussed  in  the  next  chapter.  The  procedure  to  filter  out  the 
stochastic  component  of  the  error,  described  in  this  section  does  not  apply  to  such  data 
sets  because  it  is  based  on  a measurement  model  which  assumes  a sinusoidal  steady  state. 

The  third  case  is  when  the  system  under  investigation  is  evolving  very  slowly.  In 
this  case  the  change  in  the  system  during  one  complete  scan  is  very  small  and  negligible, 
but  can  be  observed  over  the  course  of  several  experiments.  Such  a system  is  called  a 
pseudo-stationary  system  and  it  can  be  shown  that  each  scan  is  consistent  with  Kramers- 
Kronig  relations.  The  results  for  such  a pseudo-stationary  system  are  illustrated  in  Figure 
4.3.  Results  of  six  replicate  (consecutive)  experiments  for  copper  in  alkaline  solution  are 
shown  in  Figure  4.3.  Each  spectrum  shown  in  the  figure  was  found  to  be  consistent  with 
the  Kramers-Kronig  relations,  but  a careful  analysis  of  the  data  using  the  measurement 
model  approach  showed  that  the  data  were  not  replicate,  i.e.  the  system  changed  from  one 
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experiment  to  another  [25],  A measurement  model  with  eight  Voigt  elements  was 
regressed  to  the  combined  data  set  (i.e.,  a single  measurement  model  was  regressed  to  all 
six  data  sets  simultaneously).  The  normalized  residual  errors  obtained  from  the  regression 
are  shown  in  Figures  4.4  and  4.5.  Examination  of  the  real  and  imaginary  part  of  the 
residual  errors  shows  that  the  system  changed  from  one  experiment  to  the  other  and  the 
residuals  for  the  six  experiments  can  be  distinguished  from  one  another.  The  plot  of  the 
imaginary  part  (er .residual)  versus  the  real  part  (ej,reaidmi)  of  the  residuals  is  shown  in  Figure 
4.6.  The  plot  suggests  that  the  residuals  are  strongly  correlated.  Examination  of  real 
(circles)  and  imaginary  (triangles)  a (calculated  using  equation  4.3.7)  shown  in  Figure  4.7 
as  a function  of  frequency  suggests  that  they  are  not  equal.  This  observation  is  counter  to 
the  result  obtained  for  stationary  data,  where  the  real  part  of  the  standard  deviation  was 
observed  to  be  equal  to  the  imaginary  part.  Hence  for  non-stationary  systems  the 
stochastic  component  of  the  error  cannot  be  calculated  by  taking  the  standard  deviation  of 
the  residual  errors  directly.  A procedure  to  filter  out  the  non-stationary  component  of  the 
error  is  outlined  in  the  next  section.  The  significance  of  the  lines  in  Figure  4.7  will  be 
discussed  in  later  sections. 

4.2.3  Filtering  Technique 

A technique  to  filter  out  ens  from  non-stationary  data  to  obtain  estochastic  is  proposed 
in  this  section.  The  data  set  from  Figure  4.3  is  used  to  illustrate  the  technique.  For  non- 
stationary  data,  instead  of  regressing  a measurement  model  to  the  combined  data  set,  each 
data  is  regressed  separately.  The  measurement  model  parameters  for  each  data  set  are 
now  different  because  the  system  changed  from  one  experiment  to  the  other.  Hence,  by 
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regressing  a new  measurement  model  to  each  individual  data  set,  the  changes  of  the 
experimental  conditions  are  incorporated  into  the  measurement  model  parameters.  A 
consequence  s„  is  now  equal  to  zero  for  each  separate  regression.  The  normalized  real 
and  imaginary  residual  errors  for  six  regressions  are  shown  in  Figure  4.8  and  4.9.  The 
results  show  that  the  six  data  sets  cannot  be  distinguished  from  one  another.  The  plot  of 
the  imaginary  versus  real  residual  errors,  shown  in  Figure  4.10,  suggests  that  the  residuals 
are  not  correlated  anymore.  The  real  and  imaginary  standard  deviations  of  the  residual 
errors  are  shown  as  a function  of  frequency  in  Figure  4.11.  The  real  (circles)  and 
imaginary  (triangles)  ct  are  now  equal  to  each  other.  Hence  one  can  conclude  that  by 
regressing  a measurement  model  to  individual  spectra  the  non-stationary  error  s„s  can  be 
filtered  out  to  give  the  stochastic  error  ^stochastic ■ Although  the  solid  state  data  shown  in 
Figure  4.1  was  assumed  to  be  stationary  as  a first  approximation,  application  of  the 
filtering  algorithm  to  the  data  resulted  in  a better  estimation  of  s stochastic- 

4.3  Model  for  the  Error  Structure 

In  this  section  a semi-empirical  model  for  standard  deviation  of  the  error,  a is 
proposed.  The  literature  shows  that  relatively  little  work  has  been  done  in  developing  a 
proper  model  for  the  error  structure  of  impedance  data  [37-40],  Macdonald  et  al.  have 
proposed  a power-law  model  for  the  variance  v [37], 

Vr  =Ot2r  +p2r|Zr|2Yr 
v.  =ct2,+|32,|Zr|2Y' 


(4.8) 
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where  oc,P,  and  y are  parameters  of  the  model,  Zr  and  Zj  are  the  real  and  imaginary  part  of 
impedance,  vr,  Vj,  and 


vr  =ct2r 
vj=g2J 


(4.9) 


where  ar  and  ct,  are  the  real  and  imaginary  part  of  the  standard  deviation.  In  this  model  it  is 
assumed  that  the  real  and  imaginary  parts  of  the  error  are  not  equal  to  one  another. 
However  the  experimental  evidence  in  Figures  4.2  and  4. 1 1 shows  that  the  real  part  of  the 
standard  deviation  is  equal  to  the  imaginary  part. 

In  this  work,  based  on  the  experimental  observations,  it  is  assumed  that, 

CTr=a;.  =a  (4.10) 

Under  the  assumption  that  the  fundamental  impedance  measurement  in  the  instrumentation 
is  the  magnitude  | Z | and  the  phase  lag  <j>,  the  a for  the  real  and  imaginary  components  can 
be  expressed  as, 


a,  = 


a,  = 


|Z| 

dZr 

aj> 

ciz\ 

SZ, 

dZf 

ai> 

\z\ 

Jlz\ 


Jlz\ 


(4.11) 


where  o | z I &nd  are  the  standard  deviation  of  the  magnitude  and  phase  lag.  Following 
the  instrument  specifications,  as  a first  approximation,  the  error  in  the  phase  lag  was 
assumed  to  be  a constant;  whereas  the  error  in  the  magnitude  was  assumed  to  be 
proportional  to  the  signal  with  a term  added  to  account  for  the  poor  signal-to-noise  ratio 
experienced  when  there  is  mismatch  between  the  system  impedance  and  the  measuring 
resistor.  Thus, 
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a*  =a 


Izl 

ctizi=(P+yV)Iz| 


i? 


(4.12) 


where  a,  3,  and  y are  constants,  and  R„,  is  the  value  of  the  current 
measuring  resistor.  Parameters  a,  3,  and  y will  depend  on  the  specific  instruments  being 
used  for  impedance  measurements  ( i.e.,  the  values  of  these  parameters  will  change  with 
the  different  models  of  potentiostat  and  frequency  response  analyzer).  From  equation  4.11 
and  4. 12,  the  errors  in  the  real  and  imaginary  components  become 


<Jr  = ar 


Zj 


+ Pr\Zr\  + Yr~\Zr 


R 


m 


\A 


(4.13) 


<7j=aj\Zr\+PjZj  +n~\2, 

f'-m 

In  order  to  reconcile  equation  4.13  with  the  observation  that  the  real  and  imaginary 
standard  deviation  are  equal,  and  the  following  revised  error  structure  is  proposed: 

\z\2 


Cj  =cr  =o=a  |Zr|  + 3 


+ y 


R_ 


(4.14) 


Rewriting  equation  4. 14  in  polar  co-ordinates  we  obtain 


\z  I z 

°‘=(aiir+p^+TT-^l+N) 

I Z, 


Izl. 


(4.15) 


a^=(aii+^+yK'x^+z') 


Comparison  of  equation  4.15  with  equation  4.12  suggests  that  the  initial  assumptions 
about  the  errors  in  phase  lag  and  the  magnitude  are  incorrect  and  that  the  errors  in  phase 
lag  are  not  constant. 
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Values  for  a,  P,  and  y were  calculated  for  impedance  data  collected  on  a 
potentiostat  SI1286  and  a frequency  response  analyzer  FRA1250.  The  data  used  for  the 
estimation  of  the  error  structure  parameters  were  collected  for  n-GaAs  at  temperatures 
ranging  from  320K  to  420K.  The  frequency  range  used  for  the  experiment  was  1Hz  to 
65000Hz.  The  wide  range  of  temperature  and  frequency  ensured  a wide  range  of 
impedance  values.  Five  replicate  experiments  were  conducted  at  each  temperature.  The 
results  of  the  experiments  at  320K  are  shown  in  Figure  4.1.  As  a first  approximation  it  was 
assumed  that  the  system  was  stationary  and  the  stochastic  component  of  the  error  stochastic 
and  its  standard  deviation,  a were  calculated  for  data  sets  for  the  GaAs  sample  at  320, 
340,  360,  380,  and  400K,  as  described  in  section  4.2.1.  Equation  4.14  was  then  regressed 
to  the  standard  deviation  values  to  obtain  the  values  of  a,  P,  and  y.  The  parameter 
estimates  obtained  are  given  in  Table  4.1.  The  results  of  the  regression  for  GaAs  at  320K 
are  shown  in  Figure  4. 12  and  the  results  for  400K  are  presented  in  Figure  4. 13.  The  circles 
and  triangles  represent  the  real  and  imaginary  part  of  the  standard  deviation,  respectively, 
and  the  solid  line  represents  the  model  for  the  error  structure  given  by  equation  4.14.  The 
dashed  lines  represent  the  contribution  of  the  different  terms  in  equation  4. 14. 

The  filtering  algorithm  described  in  section  4.2.2  was  applied  to  the  above  data  set 
to  give  a better  estimate  for  the  stochastic  errors.  Equation  4.14  was  regressed  to  the 
filtered  errors  and  the  parameter  for  the  error  structure  model  are  shown  in  Table  4. 1 . The 
filtered  errors  for  GaAs  at  320K  and  400K,  with  the  new  model,  are  shown  in  Figure  4. 14 
and  4. 15  respectively.  The  new  standard  deviations  (Figure  4. 14)  are  now  smaller  than  the 
unfiltered  case  (Figure  4.12)  and  the  real  and  imaginary  parts  of  the  standard  deviations 
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are  closer.  The  regression  results  in  Figure  4.14  and  4.15  show  that  the  agreement  of  the 
model  for  the  error  structure  (equation  4.14)  to  the  experimental  standard  deviations  is 
quite  good.  Similar  agreement  was  found  at  all  temperatures.  The  parameters  from  Table 
4.1,  for  filtered  errors,  were  used  to  predict  the  errors  for  corrosion  of  copper  data  shown 
in  Figures  4.3.  The  results  are  presented  in  Figure  4.11.  The  model  shows  a good 
agreement  with  the  experimentally  obtained  standard  deviations.  The  validity  of  the 
equation  is  supported  by  the  fact  that  a three-parameter  model  provides  a good  agreement 
for  solid  state  systems  as  well  as  for  corrosion  systems,  for  data  collected  under  a wide 
variety  of  experimental  conditions  and  for  errors  ranging  from  milli-ohm  to  mega-ohm  in 
magnitude. 


4.4  Applications  of  the  Error  Structure 

A priori  knowledge  of  the  error  structure  for  impedance  experiments  helps  in  the 
analysis  of  impedance  data  and  the  design  of  experiments.  Identification  of  the  various 
components  of  the  error  can  help  guide  the  design  of  experiments  because  one  can  thus 
minimize  the  source  of  error.  Furthermore,  the  model  for  the  error  structure  can  also  be 
used  to  weight  data  during  regression  of  a model  to  impedance  data. 

Proper  weighting  of  data  is  crucial  during  non-linear  regression  of  a model  to 
impedance  data  to  get  unbiased  estimates  of  the  parameters.  In  the  absence  of  explicit 
weighting,  the  objective  function  J,  minimized  by  the  regression  algorithm  for  complex 
fitting,  is  given  by 


•'=S(z„-z„)J  + Z (zM-z,t)J 

* k 


(4.16) 
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where  Zr>k  and  ZJ>k  represent  the  real  and  imaginary  part  of  the  data,  respectively,  and 

A /y 

Zr  t and  ZJk  represent  the  real  and  imaginary  part  of  the  model.  The  objection  to 
approach  in  equation  4.16  is  that  the  largest  magnitude  terms  contribute  the  greatest 
amount  to  the  objective  function,  and  potentially  important  high-frequency  information 
(where  Zj  tends  towards  zero  and  Z,  becomes  small)  can  often  have  a negligible 
contribution  to  equation  4. 16,  and  hence  such  data  is  effectively  ignored. 

The  most  commonly  used  weighting  strategy  found  in  the  literature  is  proportional 
weighting,  either  using  the  experimental  data: 


j=Y. 


* z,/  r V 

or,  as  suggested  by  Macdonald  and  Potter  [37],  using  model  estimates: 


j=1L—-  z,j)' 

k Zr,k 


(4.17) 


(4.18) 


* V 

Proportional  weighting  gives  equal  weights  to  all  data  points.  The  problem  with  this 
approach  is  that,  data  points  with  “high  noise”  content  contribute  equally  to  the  regression 
as  data  points  with  “low  noise”  content. 

The  most  defensible  weighting  approach  is  to  weight  the  regression  by  the 
measured  variance  of  the  data  a2,  i.e. 


T ^(Zr,k  Zr,k ) (Zj,k  Zj,k )2 

J=L — —2 — +2-- 


Jr,k 


(4.19) 


where  crrjc  and  aJ>k:  are  the  real  and  imaginary  part  of  the  standard  deviation  at  each 
frequency.  The  use  of  variance  to  weight  data  ensures  that  data  points  with  “low  noise” 


content  are  emphasized  and  the  data  points  with  “high  noise”  content  are  de-emphasized. 
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Replicate  experiments  are  needed  to  calculate  ct2  experimentally,  however,  in  practice 
sometimes  it  is  not  possible  to  obtain  replicate  experiments.  If  the  experimental  values  for 
the  variance  are  not  available,  a model  for  the  variance  a2  could  be  used: 


r Zr,k-Zr,k)2  ^(ZM-ZM)a 

J=L — — +2- — — y 


r.k 


(4.20) 


The  model  for  the  error  structure  given  by  equation  4.14  provides  an  accurate  estimate  for 
the  variance  and  can  be  used  to  weight  the  regression. 


4.5  Conclusions 

In  this  chapter  the  different  errors  that  occur  in  impedance  measurements  were 
discussed.  Algorithms  to  obtain  the  stochastic  component  of  the  error  from  stationary  and 
non-stationary  data  were  outlined.  A semi-empirical  model  for  the  error  structure  was 
proposed.  The  error  structure  proposed  in  this  chapter  is  preliminary  and  work  needs  to 
be  done  to  further  refine  the  model  to  incorporate  other  sources  of  errors.  It  was  proposed 
that  the  model  for  the  error  provides  a proper  way  of  weighting  data  during  regression. 
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Table  4.1:  Parameter  Estimates  for  the  Error  Structure  Model  (Equation  4.14)  Obtained 
by  Regression  to  replicate  Data  for  n-GaAs. 


Parameter 

Unfiltered 

Filtered 

a 

3.27  x 10"4  ± 8.8  x 10'5 

8.12  x 10-4  ± 9.07  x 10'5 

P 

1.20  x 10'3  ± 7.5  x 10’5 

9.33  x 10^1 4.99  x 10"4 

Y 

2.83  x 10-4  ± 1.5  x 10“* 

2.31  x 10"4  ± 7.88  x 10^ 
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Figure  4.1:  Five  replicate  measurements  for  an  n-GaAs  single  crystal  held  at  320K. 
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Figure  4.2:  Unfiltered  standard  deviation  for  GaAs  single  crystal  at  320K.  Circles 
represent  the  real  and  the  triangles  represent  the  imaginary  part  of  the  standard  deviation. 
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Figure  4.3:  Six  replicate  Impedance  measurements  for  Copper  in  alkaline  solution. 
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Figure  4.4:  Real  residual  errors  for  the  regression  of  a single  measurement  model  to  the 
data  shown  in  Figure  4.3. 
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Figure  4.5:  Imaginary  residual  errors  for  the  regression  of  a single  measurement  model  to 
the  data  shown  in  Figure  4.3. 
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Figure  4.6:  Plot  of  the  imaginary  errors  versus  the  real  unfiltered  errors  for  the  data 
shown  in  Figure  4.3. 
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Figure  4.7:  Real  and  imaginary  unfiltered  errors  for  the  data  shown  in  Figure  4.3  as  a 
function  of  frequency.  The  circles  represent  the  real  and  the  triangles  represent 
the  imaginary  errors. 


92 


q 2 — rTmTOj — rmrnr| — rrrmnj — rTTnnt| — rrrnnj| — m 


M 


M 


M 


■0.2 


i i m mi  i i i mill  i i i mill  i i i mill  i j. UMsL.  LJU 


UUll 


Figure  4.8:  Real  residual  errors  for  the  regression  of  a single  measurement  model  to  the 
data  shown  in  Figure  4.3. 
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Figure  4.9:  Imaginary  residual  errors  for  the  regression  of  separately  measurement  model 
to  the  data  shown  in  Figure  4.3. 


94 


er,  n 


Figure  4.10:  Plot  of  the  imaginary  errors  versus  the  real  errors  for  the  data  shown  in 
Figure  4.3. 
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Figure  4.11:  Real  and  imaginary  errors  for  the  data  shown  in  Figure  4.3  as  a 
function  of  frequency.  The  circles  represent  the  real  and  the  triangles  represent 
the  imaginary  errors. 
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Figure  4.12:  The  solid  line  represents  the  preliminary  error  structure  for  the  n-GaAs 
sample  held  at  320  K.  The  dashed-dot  line  represents  the  contribution  of  the  imaginary 
term  to  the  error.  The  dashed-dot-dot  line  represents  the  contribution  of  the  real  term  and 
the  dashed  line  is  the  measuring  resistor  term. 
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Figure  4.13:  The  solid  line  represents  preliminary  error  structure  for  the  n-GaAs  sample 
held  at  400  K.  The  dashed-dot  line  represents  the  contribution  of  the  imaginary  term  to  the 
error.  The  dashed-dot-dot  line  represents  the  contribution  of  the  real  term  and  the  dashed 
line  is  the  measuring  resistor  term. 
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Figure  4.14:  The  solid  line  represents  modified  error  structure  for  the  n-GaAs  sample  held 
at  320  K.  The  dashed-dot  line  represents  the  contribution  of  the  imaginary  term  to  the 
error.  The  dashed-dot-dot  line  represents  the  contribution  of  the  real  term  and  the  dashed 
line  is  the  measuring  resistor  term. 
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Figure  4.15:  The  solid  line  represents  modified  error  structure  for  the  n-GaAs  sample  held 
at  400  K.  The  dashed-dot  line  represents  the  contribution  of  the  imaginary  term  to  the 
error.  The  dashed-dot-dot  line  represents  the  contribution  of  the  real  term  and  the  dashed 
line  is  the  measuring  resistor  term. 


CHAPTER  5 


MEASUREMENT  MODELS:  CONSISTENCY  WITH 
KRAMERS-KRONIG  RELATIONS 

In  principle,  the  Kramers-Kronig  relations  can  be  used  to  determine  whether  the 
impedance  spectrum  of  a given  system  has  been  influenced  by  time-dependent  phenomena. 
Although  this  information  is  critical  to  the  analysis  of  impedance  data,  the  Kramers-Kronig 
relations  have  not  found  widespread  use  in  the  analysis  and  interpretation  of 
electrochemical  impedance  spectroscopy  data  due  to  difficulties  with  their  application.  The 
integral  relations  require  data  for  frequencies  ranging  from  zero  to  infinity,  but  the 
experimental  frequency  range  is  often  constrained  by  instrumental  limitations  or  by  noise 
attributable  to  the  instability  of  the  electrode. 

The  Kramers-Kronig  relations  have  been  applied  to  electrochemical  systems  by 
direct  integration  of  the  equations,  by  experimental  observation  of  stability  and  linearity, 
and  by  regression  of  electrical  circuit  models  to  the  data.  Each  of  these  approaches  has  its 
merits  and  its  disadvantages.  In  this  chapter,  the  Kramers-Kronig  equations  and  the 
methods  used  to  apply  them  to  electrochemical  impedance  spectra  will  be  reviewed.  The 
disadvantages  associated  with  current  methods  used  to  check  experimental  data  for 
consistency  with  the  Kramers-Kronig  relations  can  be  circumvented  by  application  of 
measurement  models  to  impedance  spectra. 
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5 . 1 The  Kramers-Kronig  Relations 


The  Kramers-Kronig  relations,  developed  for  the  field  of  optics,  are  integral 
equations  which  constrain  the  real  and  imaginary  components  of  complex  quantities  for 
systems  that  satisfy  conditions  of  causality,  linearity,  and  stability  [41-44],  Bode  [45] 
extended  the  concept  to  electrical  impedance  and  tabulated  various  forms  of  the  Kramers- 
Kronig  relations.  Several  transformations  used  in  electrochemical  literature  are  given 
below.  The  imaginary  part  of  the  impedance  can  be  obtained  from  the  real  part  of  the 
impedance  spectrum  through 


where  Z^co)  and  Zj(co)  are  the  real  and  imaginary  components  of  the  impedance  as 
functions  of  frequency  co.  The  real  part  of  the  impedance  spectrum  can  be  obtained  from 
the  imaginary  part  through 


if  the  high  frequency  asymptote  for  the  real  part  of  the  impedance  is  known,  and  through 


if  the  zero  frequency  asymptote  for  the  real  part  of  the  impedance  is  known.  A relationship 
between  the  phase  angle  <J>(co)  and  modulus  | Z(co)  | is  also  available,  e.g., 


(5.1) 


(5.2) 


(5.3) 


(5.4) 
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In  response  to  the  integration  limits,  a fourth  constraint,  that  the  impedance  approach 
finite  values  at  frequency  limits  of  zero  and  infinity,  is  commonly  added.  This  constraint, 
sometimes  claimed  to  prevent  application  of  the  Kramers-Kronig  relations  to  capacitive 
systems,  is  in  fact  not  needed  because  a simple  variable  substitution  [45]  can  be  used  if  the 
imaginary  part  of  the  impedance  tends  to  infinity  according  to  l/o  as  co— ko0  [46], 

5.2  Review  of  Methods  for  Determining  Consistency 

The  usual  approach  in  interpreting  impedance  spectra  is  to  regress  a model  to  the 
data.  The  models  employed  are  typically  linear  and  assume  conditions  of  a sinusoidal 
steady  state.  It  is  important,  therefore,  that  the  impedance  response  be  characteristic  of  a 
system  that  is  causal,  linear,  and  stable.  The  condition  of  linearity  can  be  achieved  by  using 
sufficiently  small  amplitude  perturbations.  The  condition  of  stability  requires  that  the 
system  return  to  its  original  condition  when  the  perturbing  signal  is  terminated.  An 
additional  implied  constraint  of  stationary  behavior  may  be  difficult  to  achieve  in 
electrochemical  systems  such  as  corrosion  where  the  electrode  may  change  significantly 
during  the  time  required  to  collect  impedance  data.  It  is,  therefore,  of  practical  importance 
to  the  experimentalist  to  know  whether  the  data  taken  does  in  fact  satisfy  the  Kramers- 
Kronig  relations.  The  approaches  taken  to  ascertain  the  degree  of  consistency  include 
direct  integration  of  the  Kramers-Kronig  relations,  experimental  replication  of  data,  and 
regression  of  electrical  circuit  analogues  to  the  data. 
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5.2.1  Direct  Integration  of  Kramers-Kronig  Relation 

The  Kramers-Kronig  relations  provide  a unique  transformation  that  can  be  used  to 
predict  one  component  of  the  impedance  if  the  other  is  known  over  the  frequency  limits  of 
zero  to  infinity.  The  usual  way  of  using  the  Kramers-Kronig  equations,  therefore,  is  to 
calculate  the  imaginary  component  of  impedance  from  the  measured  real  component 
using,  for  example,  equation  5.1  and  to  compare  the  values  obtained  to  the  experimental 
imaginary  component.  Alternatively,  the  real  component  of  impedance  can  be  calculated 
from  the  measured  imaginary  values  using  equation  5.2  or  5.3.  The  major  difficulty  in 
applying  this  approach  is  that  the  measured  frequency  range  may  not  be  sufficient  to  allow 
integration  over  the  frequency  limits  of  zero  to  infinity.  Therefore,  discrepancies  between 
experimental  data  and  the  impedance  component  predicted  through  application  of  the 
Kramers-Kronig  relations  could  be  attributed  to  the  use  of  a frequency  domain  that  is  too 
narrow,  as  well  as  to  the  failure  to  satisfy  the  constraints  of  the  Kramers-Kronig  equations. 
The  Kramers-Kronig  relations  can,  in  principle,  be  applied  with  a suitable  extrapolation  of 
the  data  into  the  unmeasured  frequency  domain.  Several  methods  for  extrapolation  have 
appeared  in  electrochemical  literature. 

Kendig  and  Mansfeld  [47]  proposed  extrapolating  an  impedance  spectrum  into  the 
low  frequency  domain  under  the  assumption  that  the  imaginary  impedance  is  symmetric, 
thus,  the  polarization  resistance  Rp  otherwise  obtained  from 


0 


(5.5) 


104 


would  be  obtained  from 


R 


7>  ~ 


(5.6) 


where  comax  is  the  frequency  at  which  the  maximum  in  the  imaginary  impedance  is 


observed.  This  approach  is  limited  to  systems  which  can  be  modeled  by  a single  relaxation 
time  constant  [48,49],  The  limitation  to  a single  time  constant  is  severe  because  multiple 
elementary  processes  with  different  characteristic  time  constants  are  usually  observed  in 
electrochemical  impedance  spectra. 

Macdonald  and  Urquidi-Macdonald  [48]  have  presented  an  approach  based  on 
extrapolating  polynomials  fit  to  impedance  data.  The  experimental  frequency  domain  was 
divided  into  several  segments,  and  the  individual  impedance  components  Zr(©)  and  Zj(co) 
were  fitted  to  a polynomial  expression  given  by 


which  was  extrapolated  into  the  unmeasured  frequency  domain.  The  Kramers-Kronig 
equation  (e.g.,  equation  (5.1  or  5.3))  was  integrated  numerically  using  the  extrapolated 
piece-wise  polynomial  fit  for  either  the  real  or  the  imaginary  component  of  the  impedance, 
respectively  [50-52],  While  the  extrapolation  algorithm  was  applied  successfully  to  a 
variety  of  systems  (including  synthetic  impedance  data  derived  from  equivalent  electrical 
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circuits  and  experimental  systems  such  as  Ti02-coated  carbon  steel  in  aqueous  HC1/KC1 
solutions),  such  extrapolation  of  polynomials  is  unreliable  over  a broad  frequency  range. 

Haili  [53]  provided  an  alternative  approach  based  on  the  expected  asymptotic 
behavior  of  a typical  electrochemical  system.  For  extrapolation  to  co=0,  the  imaginary 
component  Zj(©)  was  assumed  to  be  proportional  to  to  as  ©->0  consistent  with  the 
behavior  of  a Randles-type  equivalent  circuit.  This  approach  would  apply  as  well  to  a 
Warburg  impedance,  which  is  also  nearly  proportional  to  © as  ©—>0.  The  real  impedance 
approaches  a constant  limit  which  is  the  sum  of  the  ohmic  solution  resistance  and  the 
polarization  resistance.  The  extrapolation  in  this  region  involves  only  one  adjustable 
parameter  whose  value  will  approach  Z^w^,,)  if  ©min  is  sufficiently  small.  At  high 
frequencies,  the  imaginary  component  was  assumed  to  be  inversely  proportional  to 
frequency  as  ©->«  and  the  real  component  was  assumed  to  approach  a constant 
equivalent  to  the  ohmic  solution  resistance  R*.  The  method  of  Haili  guarantees  well- 
behaved  extrapolation  of  the  impedance  spectrum  at  upper  and  lower  frequency  limits  with 
only  five  adjustable  parameters.  Haili's  work  confirmed  the  importance  of  extrapolating 
impedance  data  to  both  zero  and  infinite  frequency  when  evaluating  the  Kramers-Kronig 
relations. 

Esteban  and  Orazem  [54,55]  presented  an  approach  which  circumvented  the 
problems  associated  with  extrapolations  of  polynomials  and  yet  avoided  making  a priori 
assumption  of  a model  for  asymptotic  behavior.  Esteban  and  Orazem  suggested  that, 
instead  of  predicting  the  imaginary  impedance  from  the  measured  real  impedance  using 
Equation  5.1  or,  alternatively,  predicting  the  real  impedance  from  the  measured  imaginary 
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values  using  Equations  5.2  or  5.3,  both  equations  could  be  used  simultaneously  to 
calculate  the  impedance  below  the  lowest  measured  frequency  comin.  The  low-frequency 
limit  coo  is  an  adjustable  parameter  that  is  typically  three  or  four  orders  of  magnitude 
smaller  than  ow  The  calculated  impedance,  in  the  domain  cc>o<co<comin,  "forces"  the 
experimental  data  set  to  satisfy  the  Kramers-Kronig  relations  in  the  frequency  domain 
coo^o^max-  The  parameter  a>0  is  chosen  to  satisfy  the  requirements  that  the  real 
component  of  the  impedance  spectrum  attains  an  asymptotic  value  and  that  the  imaginary 
component  approaches  zero  as  ©— kdo.  Internal  consistency  between  the  impedance 
components  also  requires  that  the  calculated  functions  be  continuous  with  the 
experimental  data  at  to*.  These  requirements  cannot  simultaneously  be  satisfied  by  data 
from  systems  that  do  not  satisfy  the  constraints  of  the  Kramers-Kronig  relations;  therefore, 
discontinuities  between  experimental  and  extrapolated  values  were  attributed  to 
inconsistency  with  the  Kramers-Kronig  relations.  The  approach  described  by  Esteban  and 
Orazem  [54,55]  is  different  from  other  algorithms  presented  here  because  the  Kramers- 
Kronig  relations  themselves  were  used  to  extrapolate  data  to  frequencies  below  the  lowest 
measured  frequency.  Extrapolation  of  polynomials  or  a priori  assumption  of  a model  was 
thereby  avoided. 

While  each  of  the  algorithms  described  here  have  been  applied  to  some 
experimental  data  with  success,  any  approach  toward  extrapolation  can  be  applied  over 
only  a small  frequency  range  and  cannot  be  applied  at  all  if  the  experimental  frequency 
range  is  so  small  that  the  data  do  not  show  a maximum  in  the  imaginary  impedance.  The 
extrapolation  approach  for  evaluating  consistency  with  the  Kramers-Kronig  relations 
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cannot  be  applied,  therefore,  to  a broad  class  of  experimental  systems  for  which  the 
unmeasured  portion  of  the  impedance  spectrum  at  low  frequencies  is  not  merely  part  of  a 
tail"  but  instead  represents  a significant  portion  of  the  impedance  spectrum. 

5.2.2  Experimental  Checks  for  Consistency 

Experimental  methods  can  be  applied  to  check  whether  impedance  data  conform  to 
the  Kramers-Kronig  assumptions.  A check  for  linear  response  can  be  made  by  observing 
whether  spectra  obtained  with  different  magnitudes  of  the  forcing  function  are  replicate. 
Stationary  behavior  can  be  also  be  identified  experimentally  by  replication  of  the 
impedance  spectrum.  Spectra  are  replicate  if  the  spectra  agree  within  the  expected 
frequency-dependent  measurement  error.  If  the  experimental  frequency  range  is  sufficient, 
the  extrapolation  of  the  impedance  spectrum  to  zero  frequency  can  be  compared  to  the 
corresponding  values  obtained  from  separate  steady  state  experiments.  The  experimental 
approach  to  evaluating  consistency  with  the  Kramers-Kronig  relations  shares  constraints 
with  direct  integration  of  the  Kramers-Kronig  equations.  Because  extrapolation  is 
required,  the  comparison  of  the  DC  limit  of  impedance  spectra  to  steady  state 
measurement  is  possible  only  for  systems  for  which  a reasonably  complete  spectrum  can 
be  obtained.  Experimental  approaches  for  verifying  consistency  with  the  Kramers-Kronig 
relations  by  replication  are  further  limited  in  that,  without  an  a priori  estimate  for  the 
confidence  limits  of  the  experimental  data,  the  comparison  is  more  qualitative  than 
quantitative.  A method  is  therefore  needed  to  evaluate  the  error  structure,  or  frequency- 
dependent  confidence  interval,  for  the  data  that  would  be  obtained  in  the  absence  of  non- 
stationary behavior. 
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5.2.3  Regression  of  Circuit  Analogues 

Electrical  circuits  consisting  of  passive  and  distributed  elements  satisfy  the 
Kramers-Kronig  relations  (see,  for  example,  the  discussion  in  references  [56-58]). 
Therefore,  successful  regression  of  an  electrical  circuit  analogue  to  experimental  data 
implies  that  the  data  satisfy  the  Kramers-Kronig  relations  [37,59],  This  approach  has  the 
advantage  that  integration  over  an  infinite  frequency  domain  is  not  required,  therefore  a 
portion  of  an  incomplete  spectrum  can  be  identified  as  being  consistent  without  use  of 
extrapolation  algorithms. 

Perhaps  the  major  problem  with  the  use  of  electrical  circuit  models  to  determine 
consistency  is  that  interpretation  of  a "poor  fit"  is  ambiguous.  A poor  fit  is  not  necessarily 
the  result  of  an  inconsistency  of  the  data  with  the  Kramers-Kronig  relations.  A poor  fit 
could  also  be  attributed  to  use  of  an  inadequate  model  or  to  regression  to  a local  rather 
than  global  minimum  (caused  perhaps  by  a poor  initial  guess). 

5.3  Method  Based  on  Measurement  Model 

In  this  work,  the  measurement  models  are  used  to  check  for  the  consistency  of  the 
experimental  data  with  the  Kramers-Kronig  relations.  The  methodology  proposed  here  is 
based  on  the  procedure  of  regression  of  circuit  analogues  to  determine  consistency  of 
data,  discussed  in  the  previous  section.  The  proposed  algorithm  differs  from  previous 
works  in  that,  by  including  a sufficient  number  of  lineshapes,  a measurement  model  based 
on  Voigt  circuit  elements  can  be  regressed  to  typical  stationary  electrochemical  impedance 
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data  [ see  chapter  2],  Therefore  an  inability  to  fit  an  impedance  spectrum  can  be  attributed 
to  the  failure  of  the  data  to  conform  to  the  assumptions  of  the  Kramers-Kronig  relations 
rather  than  the  failure  of  the  model.  Another  unique  feature  of  this  work  is  that  a 
quantitative  criteria  for  identifying  the  portion  of  an  impedance  spectra  inconsistent  with 
the  Kramers-Kronig  relations  is  also  proposed. 

5.3.1  Algorithm 

The  algorithm  used  to  check  data  consistency  takes  advantage  of  the  Kramers- 
Kronig  equations  (5. 1-5.3)  in  that  the  equations  relate  the  real  part  of  the  impedance 
spectra  to  the  imaginary  part,  i.e.  the  real  part  of  the  spectra  is  constrained  by  the 
imaginary  part  and  vice  versa.  The  measurement  model  is  regressed  to  the  real  (or 
imaginary)  part  of  the  spectra  and  the  parameters  obtained  from  the  regression  are  used  to 
predict  the  imaginary(or  real)  part  of  the  spectra.  Experimental  data  always  have  random 
or  stochastic  errors  associated  with  it.  The  presence  of  these  random  errors  gives  rise  to 
an  uncertainty  in  the  prediction  of  parameters  in  a regression.  This  uncertainty  in  the 
parameter  estimation  is  quantified  by  the  standard  deviation(a)  of  the  parameters,  i.e.  one 
can  say  with  95.4%  certainty  that  the  parameter  estimates  lie  within  2a  of  the  value 
calculated  by  regression.  Due  to  this  uncertainty  in  the  parameter  estimation,  there  is 
uncertainty  in  any  prediction  that  is  made  using  these  parameters.  The  Monte-Carlo 
simulation,  described  in  the  next  section,  is  used  to  determine  the  95.4%  confidence 
interval  for  the  prediction.  Calculation  of  the  95.4%  confidence  interval  takes  into  account 
the  stochastic  component  of  measurement  error.  Hence,  one  can  say  with  95.4%  certainty 
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that  the  data  points  which  lie  outside  this  predicted  confidence  interval  are  corrupted  by 
systematic  error,  i.e.  they  are  inconsistent. 

In  principle,  regression  can  be  performed  on  the  real  or  imaginary  part  of  the 
spectra  to  obtain  the  same  information,  but  in  practice,  due  to  stochastic  errors,  the 
information  content  of  the  real  and  imaginary  part  of  the  spectra  can  differ.  Experimental 
data  obtained  for  corrosion  of  copper  in  0.5M  Cl'  solution  with  pH  of  1 1.5  [17]  are  shown 
in  figure  5.1.  The  circles  represent  the  experimental  data,  the  solid  line  represents  the 
regression  of  a measurement  models  composed  of  eight  Voigt  circuit  elements.  The 
triangles  represent  the  deconvolution  of  the  measurement  model  into  its  eight  components. 
The  data  shown  here  exhibit  a high  frequency  asymptote,  i.e.  at  the  high  frequency  the  real 
part  of  the  impedance  tends  towards  a constant  but  finite(not  equal  to  zero)  value  and  the 
imaginary  part  tends  towards  zero.  It  was  shown  in  the  previous  chapter  that  the  absolute 
values  of  the  stochastic  error  in  the  real  and  the  imaginary  part  of  impedance  are  equal. 
The  implication  of  this  result  is  that  at  high  frequency,  as  the  imaginary  part  of  the 
impedance  goes  to  zero  and  the  real  part  goes  to  a finite,  non  zero,  value  of  the  stochastic 
error  is  finite  and  non  zero,  i.e.  the  percentage  error  in  the  imaginary  part  tends  to  infinity 
at  high  frequency.  Hence,  at  high  frequency  the  imaginary  part  of  the  impedance  has  very 
little  information  content  compared  to  the  real  part  of  the  impedance  because  the 
imaginary  part  of  the  data  is  corrupted  by  stochastic  error.  Due  to  experimental  difficulties 
(such  as  the  time  required  to  collect  low  frequency  data),  enough  data  are  not  collected  to 
show  low  frequency  asymptotic  behavior.  At  the  low  frequency  end,  since  the  stochastic 
errors  in  the  real  part  are  equal  to  stochastic  errors  in  the  imaginary  part,  both  the  real  and 


Ill 


imaginary  data  are  equally  reliable.  Examination  of  Figure  5.1  shows  that,  at  low 
frequency,  the  real  part  of  the  deconvolution  of  the  measurement  model  does  not  have 
significant  frequency  dependence.  However,  a significant  change  is  observed  in  the 
imaginary  part  of  the  deconvolution.  Hence,  we  can  conclude  that  the  imaginary  part  of  an 
impedance  spectra  usually  has  higher  information  content  at  low  frequencies. 

Based  on  these  observations  about  the  nature  of  typical  impedance  data  the 
following  algorithm  is  proposed  to  check  for  consistency  of  impedance  data: 

Case  1 : At  low  frequency . the  imaginary  part  of  impedance  is  finite  and  the  real 
part  of  the  impedance  does  not  reach  an  asymptote  and  at  high  frequency  the 
imaginary  part  of  the  impedance  soes  to  zero  and  the  real  part  of  the  impp.dnnrp 
reaches  an  asymptote. 

Step  1:  Check for  consistency  at  the  high  frequency  end. 

a)  Regress  the  measurement  model  to  the  real  part  of  the  data. 

b)  Predict  the  imaginary  part  and  the  95.4%  confidence  interval  through  Monte- 
Carlo  simulation. 

c)  High  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Step  2:  Check  for  consistency  at  the  low  frequency  end. 

a)  Regress  the  measurement  model  to  the  imaginary  part  of  the  truncated  data  set. 

b)  Predict  the  real  part  and  the  95.4%  confidence  interval  through  Monte-Carlo 
simulation. 
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c)  Low  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Case  2:  At  low  and  high  frequency,  the  imaginary  part  of  impedance  tends  to  zero 
qnd  the  real  part  of  the  impedance  reaches  an  asymptote. 

Step  1:  Check  for  consistency  at  the  high  frequency  end 

a)  Regress  the  measurement  model  to  the  real  part  of  the  data. 

b)  Predict  the  imaginary  part  and  the  95.4%  confidence  interval  through  Monte- 
Carlo  simulation. 

c)  High  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Step  2:  Check  for  consistency  at  the  low  frequency  end 

a)  Regress  the  measurement  model  to  the  real  part  of  the  truncated  data  set. 

b)  Predict  the  imaginary  part  and  the  95.4%  confidence  interval  through  Monte- 
Carlo  simulation. 

c)  Low  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Case  3:  At  low  frequency,  the  imaginary  part  of  impedance  goes  to  zero  and  the 
real  part  of  the  impedance  reaches  an  asymptote  and  at  high  frequency  the 
imaginary  part  of  the  impedance  is  finite  and  the  real  part  of  the  impedance  does 
not  reach  an  asymptote. 

Step  1:  Check for  consistency  at  the  low  frequency  end 
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a)  Regress  the  measurement  model  to  the  real  part  of  the  data. 

b)  Predict  the  imaginary  part  and  the  95.4%  confidence  interval  through  Monte- 
Carlo  simulation. 

c)  Low  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Step  2:  Check  for  consistency  at  the  high  frequency  end. 

a)  Regress  the  measurement  model  to  the  imaginary  part  of  the  truncated  data  set. 

b)  Predict  the  real  part  and  the  95.4%  confidence  interval  through  Monte-Carlo 
simulation. 

c)  High  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Case  4:  At  low  and  high  frequency,  the  imaginary  part  of  impedance  is  finite  and 
the  real  part  o f the  impedance  does  not  reach  an  asymptote. 

Step  1:  Check  for  consistency  at  the  high  frequency  end 

a)  Regress  the  measurement  model  to  the  imaginary  part  of  the  data. 

b)  Predict  the  real  part  and  the  95.4%  confidence  interval  through  Monte-Carlo 
simulation. 

c)  High  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

Step  2:  Check  for  consistency  at  the  low  frequency  end 

a)  Regress  the  measurement  model  to  the  imaginary  part  of  the  truncated  data  set. 
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b)  Predict  the  real  part  and  the  95.4%  confidence  interval  through  Monte-Carlo 
simulation. 

c)  Low  frequency  data  that  lie  outside  the  95.4%  confidence  interval  are  deemed 
inconsistent.  Delete  the  inconsistent  data  points. 

The  application  of  the  algorithm  to  experimental  data  is  described  in  the 
subsequent  sections. 

5.3.2  Monte-Carlo  Simulation 

Accurate  calculation  of  the  95.4%  confidence  interval  for  a prediction  is  crucial  for 
the  success  of  the  above  algorithm.  This  calculation  is  done  through  the  use  of  Monte- 
Carlo  simulation  [60],  Proper  understanding  of  the  concepts  underlying  Monte-Carlo 
simulation  is  very  important  to  understand  the  methodology  outlined  here.  Let  us  assume 
that  atn«  is  a parameter  vector  which  describes  data  which  are  not  corrupted  by 
stochastic  errors.  Regression  of  a model  to  real  experimental  data  Zi  results  in  the 
parameter  estimate  ai  with  the  standard  deviation  vector  a[a].  The  significance  of  this 
result  is  that  one  can  say  with  95.4%  certainty  that  a^  lies  somewhere  in  the  range 
ai±2a[a].  The  problem  now  is  to  calculate  the  corresponding  95.4%  confidence  interval 
for  Z^.  To  make  this  calculation  an  assumption  has  to  be  made  regarding  the  probability 
distribution  of  the  stochastic  errors.  In  this  work,  stochastic  errors  are  assumed  to  be 
normally  distributed.  Because  of  this  assumption,  a series  of  parameters  ai  are  generated 
such  that  the  probability  distribution  for  a^-a,  follows  a normal  distribution  with  a 
standard  deviation  of  a.  For  each  ai  there  is  a simulated  Zi.  A very  large  number  of  these 
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simulations  is  performed  and  the  standard  deviation  a[Z]  of  these  simulations  are 
calculated.  One  can  now  say  with  95.4%  confidence  that  Ztn*  lies  in  the  range  Zi±2a[Z], 
The  regression  of  a measurement  model  to  experimental  data  yields  mean  value  of 
the  circuit  parameters  (Rsoi,  A,  and  ?,•  refer  to  equation  3.1)  and  the  standard  deviation 
{c[RsoiJ,  <y[AJ  and  a[Xi])  of  the  parameters.  A large  mirck>er(nrand)  of  random 
parameters  with  gaussian  distribution  of  mean(  Rsoj , A,-  and  ?,  ) and  standard  deviation 
((j[RsoiJ,  cr[Ai]  and  <j[ r,J)  were  generated  by, 

R sol,j  ~ Rsol  Ran j (0,l)cr[./?£O/  ] 

Aik  = A,-  + itow*(0,l)a[A/]  (5.9) 

Tjj  = *i  + Rani(Q,\)a[Ti] 

where  Ranj(0,l),  Rant(0,l)  and  Rani(0,l)  are  three  unique  random  deviates  with  gaussian 
distribution  with  zero  mean  and  a standard  deviation  of  one.  A value  for  the  impedance,  at 
one  frequency,  and,  for  one  Monte-Carlo  simulation,  was  calculated  by, 


Zp  (&)  - Rsol,p 


+ 


nose 

It 

/=i 


+ jTipO) 


(5.10) 


where  nose  is  the  number  of  lineshapes  in  the  measurement  model.  The  calculation  at  each 
frequency  was  repeated  nrand  times.  A large  value  of  nrand  is  desirable  for  an  accurate 
simulation.  The  mean  and  standard  deviation  of  the  impedance  were  calculated  by, 


Z(CD)  = 


nrand 

Z 


p= 1 


nrand 


(5.11) 
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Jnrand  _ n 

I (Zp(o>)  - Z(c,))2 


(5.12) 


The  95.4%  confidence  limit  is  thus  calculated  by, 


Zup(p})  = Z(eu)  + 2<j[Zp(a)\ 
Zio(Q))  = Zia)  - 2a[Zp(co)\ 


(5.13) 


where  Zup(co)  and  Zi0(o)  are  the  upper  and  lower  confidence  intervals  respectively. 


5.3.3  Application  to  Experimental  Data 

The  above  algorithm  was  applied  to  electrochemical  impedance  data  obtained  for 
corrosion  of  copper  in  0.5M  Cl'  with  pH  11.5  shown  in  figure  5.1.  Examination  of  the 
data  shows  that  the  imaginary  part  of  the  impedance  does  not  show  a maximum. 
Therefore,  as  discussed  in  the  previous  sections,  methods  based  on  direct  integration  of 
the  Kramers-Kronig  relations  cannot  be  applied  to  check  for  consistency.  The  algorithm 
developed  by  Esteban  et.  al  [54]  was  applied  to  the  data  set  and  it  indicated  that  most  of 
the  data  set  is  inconsistent.  It  will  be  shown  later  in  this  section  that  this  is  an  erroneous 
conclusion. 

Results  of  a complex  regression  of  a measurement  model  to  the  experimental  data 
are  shown  in  Figures  5.2  and  5.3.  The  fit,  shown  here,  was  obtained  by  weighting  the  data 
with  the  model  for  the  error  structure  proposed  in  the  previous  section  and  by  using  seven 
Voigt  circuit  elements.  This  is  the  best  fit  that  could  be  obtained  for  the  data.  In  Figure  5.2 
the  circles  represent  the  experimental  data  points  and  the  solid  line  represents  the 
regression  of  the  measurement  model  to  the  experimental  data.  In  Figure  5.3  the  circles 
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represent  the  fractional  residual  errors  obtained  from  the  regression.  Examination  of  the 
residual  errors  shows  that  errors  in  the  real  part  are  about  10%  and  errors  in  the  imaginary 
part  are  about  8%.  Qualitative  examination  of  the  fit  would  suggest  that  the  regression 
seems  quite  good.  Therefore,  based  on  the  traditional  approach  of  regression  of  circuit 
analogues  to  impedance  data,  one  would  conclude  that  the  data  is  consistent. 

To  check  for  consistency  of  the  data  using  the  algorithm  proposed  in  this  work, 
the  measurement  model  was  regressed  to  the  real  part  of  the  data.  The  best  fit  was 
obtained  by  using  eight  Voigt  circuits.  The  results  of  this  regression  to  the  real  part  of  the 
data  are  shown  in  Figures  5.4  and  5.5.  The  circles  in  Figure  5.4  represent  the  experimental 
data.  The  middle  line  in  Figure  5.4a  represents  the  regression  of  the  measurement  model 
to  the  real  part  of  the  data  and  the  middle  line  in  Figure  5.4b  represents  the  model 
predictions  using  parameters  obtained  from  the  regression.  The  two  outer  lines  in  figure 
5.4a  and  5.4b  represent  the  95.4%  confidence  interval  lines  calculated  by  Monte-Carlo 
simulation.  The  residual  errors  are  shown  in  Figure  5.5.  The  circles  in  Figure  5.5  represent 
the  normalized  residual  errors  and  the  solid  lines  represent  the  95.4%  confidence  interval. 
Examination  of  the  high  frequency  end  of  the  imaginary  residuals  (Figure  5.5b)  shows  that 
the  data  points  with  frequencies  greater  than  10000  Hz  he  outside  the  95.4%  confidence 
interval.  Hence  these  points  are  deemed  to  be  inconsistent  and  are  deleted. 

The  measurement  model  is  then  regressed  to  the  imaginary  part  of  the  truncated 
data  set  and  the  real  part  of  the  impedance  and  the  95.4%  confidence  intervals  are 
calculated.  The  best  fit  was  obtained  with  eight  Voigt  circuits.  The  results  of  the 
regression  are  shown  in  Figures  5.6  and  5.7.  The  examination  of  the  residual  errors  and 
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the  95.4%  confidence  interval  for  the  real  part  in  Figure  5.7a  shows  that  two  low 
frequency  data  points  lie  outside  the  confidence  interval  and  hence  the  these  points  are 
deemed  to  be  inconsistent  and  deleted.  For  the  calculation  of  the  95.4%  confidence 
interval  five  thousand  Monte-Carlo  simulations  were  performed  at  each  frequency  to 
ensure  that  impedance  values  follow  a Gaussian  distribution.  Results  of  simulations  at  a 
frequency  of  5.2Hz  are  shown  in  Figure  5.8  to  illustrate  that  the  calculated  impedance  has 
a Gaussian  distribution.  In  Figure  5.8  the  dashed  line  represents  the  histogram  of  the  five 
thousand  simulations  performed  at  5.2Hz  and  the  solid  line  represents  the  theoretical 
normal  distribution  with  the  same  mean  and  standard  deviation.  Comparison  of  the  two 
curves  suggests  that  five  thousand  simulations  give  a good  approximation  to  a theoretical 
gaussian  distribution.  Results  for  complex  regression  of  a measurement  model  to  truncated 
data  set  are  shown  in  figure  5.9  and  5.10.  The  best  fit  for  this  regression  was  obtained  by 
using  nine  Voigt  circuits.  Examination  of  the  residual  errors  shows  that  the  maximum 
error  has  now  been  reduced  to  1.5%  for  both  real  and  imaginary  part  of  the  impedance 
compared  to  8%  for  the  complete  data  set.  It  should  be  noted  only  seven  Voigt  circuits 
could  be  regressed  to  the  complete  data  set;  but,  by  deleting  the  inconsistent  data  points 
we  were  able  to  regress  nine  Voigt  circuits  to  the  data,  thereby  increasing  the  amount  of 
information  that  can  be  extracted  from  the  data. 

In  addition  to  the  examples  of  the  algorithm  to  corrosion  of  copper  data  shown  in 
this  section,  the  approach  was  used  to  check  for  consistency  of  data  for  several  other 
systems.  The  application  of  the  algorithm  to  metal  hydride  data  is  shown  in  the  next 


chapter. 
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5.4  Conclusion 

The  use  of  measurement  models  is  superior  to  the  use  of  polynomial  fitting 
because  fewer  parameters  are  needed  to  model  complex  behavior,  and  because  the 
measurement  model  satisfies  the  Kramers-Kronig  relations  implicitly.  Experimental  data 
can,  therefore,  be  checked  for  consistency  with  the  Kramers-Kronig  relations  without 
actually  integrating  the  equations  over  frequency.  The  algorithm  proposed  in  this  chapter, 
in  conjunction  with  error  structure  weighting  provides  a robust  way  to  check  for 
consistency  of  impedance  data.  The  use  of  measurement  models  does  not  require 
extrapolation  of  the  experimental  data  set;  therefore,  inaccuracies  associated  with  an 
incomplete  frequency  spectrum  are  resolved.  For  the  application  to  a preliminary  screening 
of  the  data,  the  use  of  measurement  models  is  superior  to  the  use  of  more  specific 
electrical  circuit  analogues  because  one  can  determine  whether  the  residual  errors  are  due 
to  an  inadequate  model,  to  failure  of  data  to  conform  to  the  Kramers-Kronig  assumptions, 
or  to  experimental  noise. 
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Figure  5.1:  Deconvolution  of  a measurement  model  into  its  components.  The  circles 
represent  the  experimental  data.  The  solid  line  represents  the  measurement  model  and  the 
triangles  represent  the  deconvolution. 
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Figure  5.2:  Results  of  regression  of  a measurement  model  to  impedance  data  obtained  for 
Cu  in  0.5M  Cl'  solution  of  pH  1 1.5.  The  circles  represent  the  experimental  data.  The  solid 
line  represents  the  complex  fit  to  the  data. 
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Figure  5.3:  Residual  errors  for  regression  presented  in  Figure  5.2. 
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Figure  5.4:  Results  of  regression  of  a measurement  model  to  impedance  data  obtained  for 
Cu  in  0.5M  Cl'  solution  of  pH  11.5.  The  circles  represent  the  experimental  data.  The 
middle  line  represents  the  real  fit  to  the  data.  The  outer  lines  represent  the  95.4% 
confidence  interval  for  the  prediction. 
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(Real  Fit) 


Figure  5.5:  Residual  errors  for  regression  presented  in  Figure  5.4.  The  solid  lines  represent 
the  95.4%  confidence  interval. 
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Figure  5.6:  Results  of  regression  of  a measurement  model  to  impedance  data  obtained  for 
Cu  in  0.5M  Cl  solution  of  pH  11.5.  The  circles  represent  the  experimental  data.  The 
middle  line  represents  the  imaginary  fit  to  the  data.  The  outer  lines  represent  the  95.4% 
confidence  interval  for  the  prediction. 
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Figure  5.7:  Residual  errors  for  regression  presented  in  Figure  5.6.  The  solid  lines  represent 
the  95.4%  confidence  interval. 
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Figure  5.8:  The  dashed  line  represents  the  histogram  of  the  five  thousand  simulations 
performed  at  5.2Hz  and  the  solid  fine  represents  the  theoretical  normal  distribution  with 
the  same  mean  and  standard  deviation. 
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Figure  5.9:  Results  of  regression  of  a measurement  model  to  truncated  impedance  data. 
The  circles  represent  the  experimental  data.  The  solid  line  represents  the  complex  fit  to  the 
data. 


Figure  5. 10:  Residual  errors  for  regression  presented  in  Figure  5.9. 


CHAPTER  6 


METAL  HYDRIDES:  ELECTROCHEMICAL  IMPEDANCE 
EXPERIMENTS  AND  ANALYSIS 

6.1  Background 

Metal  hydride  electrode  have  gained  considerable  popularity  as  anodes  in 
rechargeable  batteries  [61-64],  The  metal  hydride  electrodes  have  several  advantages  over 
other  negative  electrodes.  In  these  other  systems,  the  negative  electrode  (Cd,  Zn,  Pb,  Fe 
etc.)  is  typically  made  from  relatively  pure  elemental  metals.  The  charge  and  discharge 
cycles  in  these  electrodes  is  carried  by  the  dissolution-precipitation  in  the  aqueous 
electrolyte  solutions.  The  dissolution-precipitation  mechanism  involves  the  intermediate 
formation  of  a dissolved  species,  i.e.  metal  ion  or  metal  ion  complex,  which  precipitates 
from  the  solution  to  form  a new  solid  phase.  The  electrolyte  for  these  cells  not  only  serves 
as  an  ionic  conductor,  but  also  participates  in  the  electrode  reactions.  For  example  the 
charge  process  of  a cadmium  electrode  in  alkaline  solution  is  given  by: 

Cd(OH)2  +20H'  -» [Cd(OH)4  ]2~ 

Cd(OH)J2-  +2e  ->Cdl  +40H'  (61) 

and  the  discharge  process  for  the  lead  dioxide  electrode  is  given  by: 

Pb02  +4H+  +2e'  -»  Pb2+  +2H20 

Pb2+  + HSO-4  ->  PbS04  i +H+  (6'2) 
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The  problems  with  these  electrodes  can  thus  be  summarized  as  follows  [65-67]: 

• nucleation  and  growth  of  a new  solid  phase  at  the  surface  results  in  the 
passivation  of  the  active  surface  by  a thin  layer  of  dense  fine-grained  discharged 
material. 

• recrystallization,  especially  after  long  waiting  periods,  leads  to  the  formation  of 
larger  grains  or  stable  modifications  to  the  surface  which  result  in  memory  effects. 

• dendrite  formation  due  to  mass  transport  limitations  causes  failure  by  internal 
short  circuiting. 

• participation  of  the  electrolyte  in  the  electrode  reaction  results  in  the 
conductivity  of  the  solution  being  a function  of  the  state  of  charge  of  the 
electrode. 

The  metal-hydride  electrode,  by  contrast,  uses  a chemical  reaction  that  reversibly 
incorporates  hydrogen  into  a metal  alloy.  The  negative  electrode  allows  electrochemical 
storage  and  release  of  hydrogen  during  the  charge  and  discharge  cycle  of  the  battery.  A 
nickel  hydroxide  electrode  is  used  as  the  positive  electrode  in  the  battery.  The  half  cell 
reaction  for  the  charge  and  discharge  cycles  can  be  written  as  follows: 

M + H20  + e <=>  MH  + OH 
Ni(OH)2  + OH'  <=>  NiOOH  + H20  + e' 

where  M is  the  metal  hydride  electrode.  From  equation  6.3  we  can  see  that  there  is  no  net 
change  in  the  electrolyte  concentration  over  the  charge-discharge  cycle.  This  result 
contrasts  with  other  systems,  and  the  consequence  of  this  result  is  that  the  overall  stability 
of  the  electrode  is  much  better  than  other  electrodes.  Also,  the  chemical  states  of  the 
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electrode  are  metallic,  for  both  oxidation-reduction  reaction,  so  the  electrical  conductivity 
is  high  in  both  the  charged  and  discharged  states.  The  hydrogen  atom  enters  the  metal 
lattice  during  the  formation  of  the  hydride  (charging  of  the  electrode).  Due  to  the  small 
size  of  the  hydrogen  atom  the  volumetric  expansion  of  the  lattice  is  about  10%  and  there 
are  no  crystallographic  changes  as  are  seen  in  the  oxidation-reduction  of  Cd,  Zn,  Pb,  or  Fe 
electrodes. 

Typically,  there  are  two  types  of  metal  hydride  electrodes  that  are  used 
commonly:  the  titanium-based  AB2  type  electrode  (e.g.,  TiNi2),  and  lanthanum-based  AB5 
type  electrode  (LaNis).  In  this  chapter  only  the  AB5  type  of  electrodes  will  be  discussed. 
The  metal  hydride  materials  used  for  the  battery  need  to  satisfy  an  extensive  list  of 
requirements.  These  properties  are  obtained  by  considering  a range  of  alloys  for  the 
electrode  materials  containing  elements  that,  if  used  alone,  would  be  unacceptable  for 
several  reasons.  The  elements  that  are  used  to  design  a metal  hydride  electrodes  are  Li,  C, 
Mg,  Al,  Si,  Ca,  Ti,  V,  Cr,  Mn,  Fe,  Co,  Ni,  Cu,  Y,  Zr,  Nb,  Mo,  Sn,  La,  W,  and  Re. 
Different  metals  in  this  list  affect  different  properties  of  the  electrode. 

The  most  important  property,  the  amount  of  hydrogen  the  material  can  absorb, 
determines  the  electrochemical  storage  capacity  and  hence  the  energy  storage  capacity  of 
the  battery.  The  hydrogen  storage  capacity  of  the  alloy  is  increased  by  increasing  the 
compositional  and  structural  disorder  in  the  electrode  and  thus  providing  lattice  sites  for 
the  incorporation  of  the  hydrogen  atom.  The  elements  Mg,  Ti,  V,  Zr,  Nb,  and  La  provide 
the  structural  and  compositional  disorder  and  alloys  can  be  designed  which  can  store 
hydrogen  up  to  2.5%  of  the  weight  of  the  alloy. 
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The  rechargeable  metal  hydride  electrode  operates  on  the  principle  of  the 
reversible  absorption  of  hydrogen  into  the  electrode.  To  ensure  reversibility,  the  metal  to 
hydrogen  bond  strength  is  an  important  design  parameter.  If  the  metal  hydrogen  bond 
strength  is  too  weak,  hydrogen  will  not  react  with  the  alloy  and  the  alloy  will  evolve  gas. 
If  the  bond  strength  is  too  large,  the  metal  hydride  electrode  is  extensively  oxidized  and 
does  not  store  hydrogen  reversibly.  The  bond  strength  of  the  alloy  is  adjusted  by  using 
elements  of  different  bond  strengths.  The  bond  strength  of  elemental  La,  Ti,  Zr  and  V is 
high  and  Ni  bonds  weakly  with  hydrogen. 

The  metal  hydride  batteries  operate  in  a strongly  oxidizing  medium  composed  of 
high-concentration  alkaline  electrolyte.  In  this  aggressively  oxidizing  environment, 
oxidation  and  corrosion  resistance  of  the  metal  hydride  electrode  is  critical.  Lanthanum,  Ti 
and  Zr  form  thick  dense  passive  oxides  in  alkaline  solutions  whereas  Ni  is  resistant  to 
oxidation.  The  combination  of  these  elements  makes  the  electrode  resistant  to  oxidation, 
but,  at  the  same  time,  the  oxide  films  that  are  formed  on  the  surface  contain  regions  of 
metallic  Ni.  These  regions  help  provide  the  necessary  electrical  conductivity  and  catalytic 
activity  in  the  oxide  film. 

Several  techniques  have  been  used  in  the  literature  to  evaluate  the  several 
properties  of  the  metal  hydride  electrode  after  a new  alloy  is  designed.  The  most  common 
electrochemical  technique  is  the  characterization  of  the  electrode  by  half  or  full  cell  cycling 
[68-71],  The  problem  with  this  technique  is  that  it  is  time  consuming  and  expensive.  Also, 
half  cell  cycling  yields  capacity  decay  data  and  does  not  identify  the  electrochemical 
reactions  leading  to  failure  of  the  electrode.  The  crystal  structure  of  the  electrode  has  been 
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studied  by  numerous  researchers  by  means  of  neutron  diffraction.  The  diffusion  of 
hydrogen  in  the  electrode  has  been  studied  by  techniques  such  as  NMR  and  quasi-elastic 
neutron  scattering  (QNS).  Oxidation  of  the  metal  hydride  surface  has  been  studied  by 
LEED,  Auger  electron  spectroscopy  and  photoelectron  spectroscopy  [72-74],  Relatively 
little  work  has  been  done  to  characterize  the  metal  hydride  using  electrochemical 
impedance  spectroscopy  [75],  Electrochemical  spectroscopy  is  a non-destructive  in-situ 
technique  which  has  been  used  extensively  to  study  several  class  of  electrochemical 
systems.  Electrochemical  impedance  espectroscopy  in  conjunction  with  other  techniques 
can  be  very  powerful  in  characterizing  the  metal  hydride  electrodes,  especially  for  the 
characterization  of  corrosion,  heterogeneous  kinetics  and  diffusion  of  hydrogen  in  the 
metal  hydride  electrode. 

In  this  work  a preliminary  investigation  of  two  kinds  of  metal  hydride  electrodes 
(LaNi5  and  MmNis)  was  done  using  electrochemical  impedance  spectroscopy.  It  will  be 
shown  that  electrochemical  impedance  spectroscopy  in  conjunction  with  the  proper 
analysis  of  the  measurements  using  the  measurement  models  is  a viable  technique  for  the 
investigation  of  these  class  of  electrodes.  A preliminary  model,  which  is  based  on  the 
solution  of  the  fundamental  governing  equation  for  the  electrode  will  be  shown. 

6.2  Electrochemical  Impedance  Experiments 

Electrochemical  impedance  experiments  were  performed  on  LaNi5  and  MmNi5 
ingot  electrodes.  A PAR  273 A potentiostat  was  used  to  control  the  potential  and  measure 
the  current  and  FRA  1250  frequency  response  analyzer  was  used  to  apply  the  sinusoidal 
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perturbation.  The  instruments  were  interfaced  to  a 33Mz  386  personal  computer.  The 
software  described  in  the  chapter  2 was  used  to  control  the  experiments. 

A three-cell  configuration  was  used  to  conduct  the  experiment.  A schematic  of  the 
experimental  setup  is  shown  in  Figure  6.1.  A Hg/HgO  reference  electrode  was  used  to 
monitor  the  potential.  A large  surface  area  platinum  counter  electrode  was  used  to  allow 
for  the  passage  of  current.  The  working  electrode  was  prepared  by  soldering  a nickel  wire 
was  to  a small  piece  of  the  ingot  electrode.  The  nickel  wire  provide  the  electrical  contact 
for  the  working  electrode.  The  electrode  was  then  allowed  to  set  in  epoxy  resin  for  24 
hours.  The  surface  of  the  electrode  was  polished  to  expose  0.28  cm2  of  surface  area.  The 
reference,  counter  and  working  electrode  were  immersed  in  31  wt%  KOH  electrolyte 
solution.  The  solution  was  deareated  with  Argon  gas  throughout  the  experiment  to 
prevent  the  oxidation  of  the  surface  of  the  working  electrode.  The  open  circuit  potential  of 
the  working  electrode  was  monitored  before  the  start  of  the  impedance  experiment,  and 
impedance  scans  were  performed  when  the  open  circuit  potential  reached  a stable  steady 
state  value.  Ten  consecutive  impedance  scans  were  performed  at  several  different 
operating  conditions.  Multiple  scans  were  taken  at  each  operating  condition  to  facilitate 
the  characterization  of  measurements  with  the  measurement  model  approach  described  in 
previous  chapters.  The  open  circuit  potential  and  the  polarization  current  were  recorded 
before  and  after  each  scan  to  monitor  the  drift  in  the  open  circuit  potential  during  the 


course  of  an  impedance  scan. 
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6.2.1  LaNi5 

The  steady  state  open  circuit  potential  for  the  LaNis  electrode  was  measured  to  be 
-0.66  V vs.  Hg/HgO.  Ten  consecutive  impedance  experiments  were  performed  under 
potentiostatic  regulation  at  open  circuit  potential.  A lOmV  potential  perturbation  was 
used  for  each  experiment.  A frequency  range  of  20kHz  to  0.1Hz  was  used.  The  software 
automatically  adjusted  the  current  measuring  resistor  during  the  course  of  the  experiment. 
The  experimental  results  are  shown  in  Figure  6.2.  The  open  circuit  potential  before  and 
after  each  scan  is  shown  as  a function  of  the  run  number  in  Figure  6.3.  The  open  circuit 
potential  remained  relatively  constant  during  the  course  of  the  experiment.  Experimental 
data  for  ten  impedance  experiments  under  galvanostatic  regulation,  at  zero  current,  are 
shown  in  Figure  6.4.  A target  potential  of  lOmV  was  used  as  the  perturbation  amplitude. 
The  amplitude  of  the  current  perturbation  was  adjusted  automatically  to  be  under  lOmV. 
The  open  circuit  potential  before  and  after  each  scan  is  shown  as  a function  of  the  run 
number  in  Figure  6.5.  Examination  of  Figures  6.2  and  6.4  shows  that  data  does  not  show  a 
low  frequency  asymptote.  The  absolute  value  of  the  impedance  ranges  from  about  2 Cl  at 
high  frequency  (20kHz)  to  about  10kf2  at  low  frequency  (0. 1Hz).  A set  of  ten  impedance 
experiments  were  performed  at  the  cathodic  potential  of  -0.8  vs.  Hg/HgO.  The 
experimental  data  is  shown  in  Figure  6.6.  The  frequency  range  used  for  these  experiments 
was  20kHz  to  1Hz.  The  data  was  collected  only  up  to  1Hz  because  at  cathodic  potential 
there  is  an  uptake  of  hydrogen  into  the  electrode.  This  results  in  the  data  being  non- 
stationary, especially  at  low  frequencies.  The  data  showed  both  high  and  low  frequency 
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asymptote.  The  absolute  value  of  the  impedance  ranged  from  about  2D  at  high  frequency 
(20kHz)  to  about  20D  at  low  frequency  (1Hz). 

6.2.2  MmNi5 

The  open  circuit  potential  for  the  MmNis  electrode  reached  a steady  state  at  the 
value  of -0.69V  vs.  Hg/HgO  reference  electrode.  Impedance  experiments  were  conducted 
under  potentiostatic  regulation  at  open  circuit  potential.  The  results  of  the  ten  replicate 
experiments  are  shown  in  Figure  6.7.  The  open  circuit  potential  before  and  after  each  scan 
is  shown  in  Figure  6.8.  Experiments  were  also  conducted  at  a cathodic  potential  of -1.0V 
vs.  Hg/HgO  and  at  an  anodic  potential  of  -0.5  V.  The  results  of  these  experiments  are 
shown  in  Figure  6.9  and  6.10.  The  examination  of  the  data  shows  the  presence  of  an 
inductive  loop  at  the  high  frequency  end.  It  will  be  shown  in  the  next  section  that  the 
inductive  part  of  the  spectra  is  inconsistent  with  the  Kramers-Kronig  relations  and  is 
deleted  for  the  further  analysis  of  the  data.  For  the  MmNis  the  high  frequency  impedance 
is  about  1KD.  At  low  frequency  (0.1Hz)  the  absolute  value  of  the  impedance  is  about 
1MD  for  open  circuit  experiments.  The  low  frequency  impedance  drops  to  about  0.4MQ 
for  experiments  at  cathodic  potential  and  to  about  0.6MD  for  experiments  conducted  at 
anodic  potentials. 

6.3  Analysis 

The  analysis  of  the  impedance  data  was  done  in  two  steps.  An  intermediate 
analysis  of  the  characteristics  of  the  measurement  was  done  with  the  measurement  model 
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approach.  The  information  obtained  from  the  analysis  was  used  in  the  regression  of  the 
process  model  to  data. 

6.3.1  Measurement  Model 

The  impedance  data  discussed  in  the  previous  section  was  analyzed  using  the 
measurement  models  as  discussed  in  the  previous  chapter.  The  data  was  checked  for 
consistency  with  the  Kramers-Kronig  relations  and  the  stochastic  component  of  the  errors 
in  the  experiment  was  determined. 

6.3. 1.1  LaNis 

The  impedance  data  for  LaNis  shown  in  Figure  6.2  seem  replicate.  A measurement 
model  with  eight  Voigt  circuit  elements  was  regressed  to  all  the  ten  data  sets 
simultaneously.  The  result  of  the  regression  are  shown  in  Figure  6.11.  The  circles 
represent  the  data  and  the  solid  line  represents  the  regression  of  the  measurement  model  to 
the  data.  The  residual  errors  obtained  from  the  regression  are  shown  in  Figure  6.12. 
Examination  of  the  residual  errors  shows  that  the  ten  experiments  can  be  distinguished 
from  each  other  and,  therefore,  the  data  are  not  truly  replicate.  Similar  analysis  was  done 
for  the  data  obtained  and  galvanostatic  regulation  and  similar  results  were  obtained.  After 
having  determined  that  the  system  changed  from  one  experiment  to  another  the  algorithm 
described  in  chapter  5 to  check  for  consistency  of  the  data  with  the  Kramers-Kronig 
relations  was  applied  to  the  individual  data  sets. 

The  results  of  the  regression  of  the  measurement  model  to  the  real  part  of  the  data 
and  prediction  of  the  imaginary  part  by  Monte-Carlo  simulation  are  shown  in  Figure  6.13. 
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In  Figure  6.13,  the  circles  represent  the  data  and  the  middle  line  represents  the  result  of 
the  regression.  The  outer  line  represent  the  95.4%  confidence  interval.  The  residual  errors 
as  a function  of  frequency  are  plotted  in  Figure  6.14.  Examination  of  the  high  frequency 
end  for  the  imaginary  part  of  the  residual  errors  shows  that  the  high  frequency  points 
above  10kHz  lie  outside  the  95.4%  confidence  interval  and  are  therefore  inconsistent.  The 
inconsistent  points  were  deleted  and  the  measurement  model  was  then  regressed  to  the 
imaginary  part  of  the  truncated  data  set.  The  results  of  the  regression  and  the  Monte-Carlo 
simulation  are  shown  in  Figure  6.15  and  6.16.  Examination  of  the  real  part  of  the  residual 
errors  shows  that  two  low  frequency  points  lie  outside  the  95.4%  confidence  interval  lines. 
Therefore  it  was  concluded  that  all  but  last  two  low  frequency  points  are  consistent  with 
the  Kramers-Kronig  relations.  Similar  analysis  was  done  to  all  the  ten  data  sets  obtained 
under  potentiostatic  regulation.  And  it  was  determined  that  all  the  data  points  in  the 
frequency  range  of  10kHz  and  2Hz  were  consistent  with  the  Kramers-Kronig  relations 
and  could  be  used  for  further  analysis.  Similar  analysis  was  done  to  data  obtained  under 
galvanostatic  regulation.  The  results  of  the  analysis  are  shown  in  Figures  6.17,  6.18,  6.19 
and  6.20.  For  these  data  sets  also  it  was  concluded  that  the  data  below  10kHz  was 
consistent  and  could  be  used  for  further  analysis. 

The  frequency  range  for  the  experiments  conducted  on  the  LaNi5  electrode  held  at 
a cathodic  potential  of  -0.8V  vs.  Hg/HgO  reference  electrode  was  20kHz-  1Hz.  A smaller 
frequency  range  was  used  for  these  experiments  because  the  preliminary  experiments  had 
shown  that  at  cathodic  potentials  the  system  is  non-stationary  due  to  the  uptake  of 
hydrogen  by  the  electrode.  The  data  for  the  experiments  at  this  cathodic  potential,  shown 
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m Figure  6.6,  was  checked  for  consistency  with  Kramers-Kronig  relations.  The  results  of 
the  analysis  of  one  of  the  data  sets  is  shown  in  Figures  6.21,  6.22,  6.23  and  6.24. 
Examinations  of  the  results  of  the  analysis  indicates  that  the  data  in  the  frequency  range  of 
lHz-2kHz  is  consistent.  The  high  frequency  points  above  2kHz  were  found  to  be 
inconsistent. 

The  next  step  in  the  analysis  of  the  data  using  the  measurement  model  was  the 
determination  of  the  stochastic  component  of  the  error.  As  it  was  shown  that  the 
experiments  were  not  truly  replicate,  the  stochastic  errors  cannot  be  obtained  by  taking  the 
standard  deviation  of  the  measurement,  and  the  data  had  to  be  filtered  to  determine  the 
stochastic  component  of  the  error.  The  standard  deviation  of  the  data  without  any  filtering 
was  calculated  for  the  open  circuit  experiment,  the  galvanostatic  experiments  and  the 
experiments  conducted  at  the  cathodic  potential.  The  results  are  shown  in  Figures  6.25, 
6.26  and  6.27.  The  examination  of  the  results  show  that  the  basic  assumption  about  the 
stochastic  errors,  that  the  real  part  of  the  standard  deviation  is  equal  to  the  imaginary  part 
of  the  standard  deviation,  does  not  hold.  The  stochastic  errors  for  these  experiments  were 
then  calculated  by  using  the  filtering  algorithm,  described  in  chapter  5.  The  results  of  the 
calculations  are  shown  in  Figures  6.28,  6.29  and  6.30.  These  results  seem  to  follow  the 
assumption  that  the  real  part  of  the  standard  deviation  is  equal  to  the  imaginary  part. 

6.3. 1.2  MmNij 

Examination  of  the  MmNi5  data  shown  in  Figures  6.7,  6.9  and  6.10  shows  the 
presence  of  an  inductive  loop  at  the  high  frequency  end.  The  data  was  checked  for 
consistency  with  the  Kramers-Kronig  relations  to  check  if  the  presence  of  the  inductive 
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loop  was  due  to  instrumental  artifacts  or  if  represented  a physical  phenomenon.  The 
results  of  the  analysis  of  one  data  set  collected  at  open  circuit  potential  are  shown  in 
Figures  6.31,  6.32,  6.33  and  6.34.  Examination  of  the  imaginary  residual  errors  in  Figure 
32  indicate  the  data  points  at  frequencies  higher  than  700Hz  are  inconsistent.  Whereas  the 
low  frequency  data  were  found  to  be  consistent  for  the  frequency  range  used  for  these 
experiments.  Similar  results  were  obtained  for  the  data  collected  at  anodic  and  cathodic 
potentials.  The  inconsistent  data  points  were  deleted  for  further  analysis. 

The  results  for  the  calculation  of  the  unfiltered  standard  deviation  of  the  data  are  shown  in 
Figures  6.35,  6.36  and  6.37  for  the  three  experimental  conditions.  The  results  are  similar 
to  the  results  obtained  for  the  LaNis  electrode  and  other  electrochemical  experiments. 
These  results  suggest  the  need  to  filter  the  data  to  calculate  the  true  standard  deviations. 
The  results  obtained  for  the  MmNis  electrode  after  the  application  of  the  filtering 
algorithm  are  shown  in  Figures  6.38,  6.39  and  6.40.  Examination  these  results  indicates 
the  improvement  in  the  prediction  of  the  true  standard  deviations  of  the  experiments. 

6.3.2  Process  Model 

A transport  based  mathematical  model  was  developed  for  the  system  which 
accounts  for  electrochemical  decomposition  of  water  to  form  hydrogen,  reversible 
absorption  of  hydrogen  into  the  electrode,  and  diffusion  of  hydrogen  into  the  bulk.  The 
model  was  based  on  the  following  mechanism  for  the  decomposition  of  water  at  the 
surface  of  the  electrode  to  form  hydrogen  and  hydroxide  and  the  incorporation  of 
hydrogen  into  the  bulk  of  the  metal  hydride  electrode  [61], 
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Step  1:  Electrochemical  decomposition  of  water  and  the  adsorption  of  hydrogen  at 
the  surface  of  the  electrode. 


where  M is  the  electrode  and  H is  the  hydrogen  adsorbed  on  the  surface  of  the 
electrode.  k\  is  the  rate  of  the  faradaic  reaction. 

Step  2:  Reversible  absorption  and  diffusion  of  hydrogen  into  the  bulk  of  the 
electrode. 


where  Habs  is  the  absorbed  hydrogen  and  k2  and  k3  are  the  reaction  rates  for  the 
homogeneous  absorption  and  desorption. 

Step  3 : Formation  of  the  metal  hydride. 


where  MHX  is  the  metal  hydride  phase,  x is  the  number  of  hydrogen  atoms  in  the 
metal  hydride  phase  per  molecule  of  the  base  alloy  and  k4  and  k5  are  the  rate 
constants  for  the  reversible  hydriding  reaction. 

The  model  developed  in  this  work  considered  only  the  first  two  steps  of  the  above 
mechanism.  The  governing  equation  for  the  surface  coverage  of  the  electrode  by  hydrogen 
is  given  by 


M + H20  + e~  — }M-  Hads  + OH~ 


(6.4) 


k2 

^ ~ H ads  vv  M + HafjS 

k3 


(6.5) 


k4 

M + xHafjS  s MHX 

k5 


(6.6) 


d6 

P ~^t  = k\Q  ~ ~ Vmax^1  ~ CH>}  + k3cmax(l  ~ °^CH 


(6.7) 
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where  j. 3 is  the  fraction  of  the  surface  available  for  reaction,  6 is  the  dimensionless  surface 
coverage,  Ch  is  the  hydrogen  concentration  in  the  bulk.  Cm**  is  the  maximum  hydrogen 
concentration  in  the  bulk  and  ki,  kj,  and  kj,  are  the  rate  constants  for  the  faradaic, 
absorption  and  desorption  reactions.  Tafel  kinetics  were  assumed  for  the  faradaic  reaction, 
and,  therefore,  the  rate  constant  ki  can  be  written  as 

kY  = k°iebE  (68) 

where  k i is  the  pre-exponential  factor  b is  the  Tafel  slope  and  E is  the  overpotential.  The 
faradaic  current  If  was  given  by 

IF  = FAky{\-G)  (69) 

where  F is  the  Faradays  constant  and  A is  the  surface  area  of  the  electrode. 

In  an  impedance  experiment  under  potentiostatic  regulation,  a sinusoidal 
perturbation  in  potential  (A E)  is  applied  to  the  system.  The  potential  perturbation  is  given 
by 

A E = EreJ0}t  (6.10) 

where  £ris  the  amplitude  of  the  potential  perturbation  and  a is  the  frequency.  The 

sinusoidal  perturbation  in  the  potential  results  in  a sinusoidal  response  in  the  current.  The 

current  is  a function  of  the  surface  coverage  0 (equation  6.9)  which  in  turn  is  a function  of 

the  surface  concentration  of  hydrogen  cH. 

The  hydrogen  concentration  can  be  expressed  by  the  sum  of  steady  state  and 
sinusoidal  steady  state  components, 

CH  = *H+<?H,r  +J°H,j>>W‘ 


(6.11) 
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where  c^,  r>£jj  j are  the  steady  state,  sinusoidal  real  and  imaginary  part 


respectively  of  the  total  concentration,  cH.  The  surface  coverage  can  also  be  expressed  by 
the  sum  of  the  steady  state  and  sinusoidal  steady  state  components, 


jo)t 


9 H ~ 6H  + (9H,  r + j9H, 


The  governing  equation  for  the  diffusion  of  hydrogen  in  the  metal  is  given  by 


3 CH  _ rv  ^'CH 

~ TT 


dt 


dzi 


(6.12) 


(6.13) 


The  boundary  conditions  are  given  by 


d c 


H 


D, 


d z 
d c 


= 0 


H 


d z 


H _ 


f(CH’  &) 


@Z 

@Z 


0 


and  the  initial  condition  is  given  by 

CH  ~ 0 @t  - 0 

where 


(6.14) 


(6.15) 


/(%,  0)  = *2^««(1  ~ ch)  - VmaxM  1 - 0)  (6.16) 

and  where  Dh  is  the  diffusion  coefficient  for  hydrogen  in  the  electrode.  Solutions  for  the 
steady  state  concentration  and  sinusoidal  steady  state  concentration  were  obtained 
sequentially. 


6.3.2. 1 Pseudo-Steady  state 

Under  the  assumption  of  a pseudo-steady  state,  equation  6.7  reduces  to 
q=  h+hcH(Q) 

h+k2  +*3ctf(°) 


(6.17) 
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where  c^O)  is  the  concentration  of  hydrogen  at  the  surface  of  the  electrode.  The  boundary 
condition  for  the  diffusion  equation  at  the  surface  of  the  electrode  was  modified  to 


The  unique  feature  of  the  steady  state  solution  is  the  iterative  treatment  of  the  non- 
linear boundary  conditional  the  surface  of  the  electrode  given  by  equation  6.16.  The 
boundary  condition  is  linearized  around  an  assumed  surface  concentration.  The  equations 
could  then  be  solved  analytically  to  give  a pseudo-state  concentration  profile.  The 
analytical  solution  for  the  concentration  of  hydrogen  with  linearized  boundary  condition  is 
given  by 


(6.18) 


CH  = 4^  “ HEn  exP 


(6.19) 


where 


A _ *1*2*3 

DH(kl+k1) 


(6.20) 


DH(kl+k  l) 


and 


(6.21) 


where  |tn  are  eigenvalues  evaluated  from  the  roots  of  the  equation 


Vton-r ==  = bJd^ 

>1dH 


(6.22) 
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The  resulting  concentration  profile  is  then  used  to  determine  the  a new  surface 
concentration  which  served  as  a guess  until  convergence  is  achieved..  Convergence  was 
achieved  within  five  or  six  iterations.  This  approach  failed  when  additional  processes  such 
as  nucleation  of  hydrogen  gas  bubbles  and  corrosion  were  incorporated  into  the  equations. 


6.3.2  2 Sinusoidal-steady  state 


The  impedance  response  Z of  the  system  under  sinusoidal  perturbation  is  given  by 


Z 


A E 
M 


Er  <6’23) 

~ FAkl\_{\-~e)bEr-Cer^fei)\ 

The  values  for  the  perturbations  in  the  surface  coverage  and  thus  the  perturbation  in  the 
hydrogen  concentration  at  the  surface  is  obtained  by  inserting  equation  6.11  into  equation 
6. 13  and  6. 14.  The  substitution  results  in 


<P"  Ctj  ; 

DH—Y~  = COCHyr 

oz 


D 


dlcH,r 


the  boundary  conditions  are  given  by 

@ z = 0 
-Dh 


d?H,r 


-L>h 


PcH,i-qcH>r=r 


(6.24) 


(6.25) 


where  p,  q,  r,  s,  t and  u are  functions  of  kh  k2,  h and  c«( 0),  and 
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@ z = 1 


dtHjr 

dz 

dz 


(6.26) 


The  sinusoidal-steady  state  equations  were  written  in  finite  difference  form  and 
solved  numerically  using  the  “BAND”  algorithm  developed  by  Newman  [76]  to  give  the 
values  for  the  surface  concentration.  The  value  of  the  surface  concentration  was  then  used 
to  calculate  the  faradaic  impedance.  A double  layer  capacity  was  added  in  parallel  to  the 
faradaic  impedance,  and  a solution  resistance  was  added  in  series. 


6. 3. 2. 3 Regression  of  process  model  to  data 

The  mathematical  model  was  regressed  to  the  experimental  data.  The  results  of  the 
regression  to  experimental  data  for  LaNis  electrode  at  open  circuit  under  potentiostatic 
regulation  are  shown  in  Figures  6.41  and  6.42.  The  regression  was  done  with  no 
weighting.  The  circles  in  figure  6.41  represent  the  data  and  the  solid  line  represents  the 
regression  of  the  model.  The  normalized  residual  errors  are  shown  in  Figure  6.42. 
Examination  of  Figure  6.41  shows  that  the  model  captured  the  features  of  the  data.  The 
residual  errors  in  Figure  6.42  indicate  that  while  a good  fit  is  obtained  at  the  low  and  mid 
frequency  range,  the  fit  is  poor  at  the  high  frequency  range.  The  confidence  in  the 
parameter  estimates  was  very  poor.  As  the  regression  was  done  with  no  weighting  the 
high  frequency  impedance  data  points,  which  were  small  in  magnitude  compared  to  the 
mid  and  low  frequency  points,  were  de-emphasized  in  the  regression.  Regression  of  the 
model  to  the  same  data  set  but  weighted  with  the  error  structure  proposed  in  this  work  is 
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shown  in  Figures  6.43  and  6.44.  The  circles  in  Figure  6.43  represent  the  data  and  the  solid 
line  represents  the  model.  The  normalized  residual  errors  for  the  real  and  imaginary  part  of 
the  impedance  are  shown  in  Figure  6.44.  Examination  of  the  results  shows  that  when 
proper  weighting  is  used  the  data  points  in  the  mid  and  low  frequency  range  do  not  show 
very  good  agreement  although  its  better  than  the  fit  at  the  high  frequency  end.  The  lack  of 
fit  indicates  that  additional  physical  processes,  such  as  corrosion  reaction  and  hydrogen 
evolution  reactions  and  potential  and  concentration  dependence  of  the  diffusion 
coefficient,  must  be  included  in  the  model. 

6.4  Conclusions 

In  this  chapter  it  was  shown  that  impedance  spectroscopy  is  an  effective  technique 
for  characterizing  metal  hydride  electrodes.  The  concept  of  measurement  model  was 
applied  to  characterize  the  measurements.  A process  model  based  on  the  numerical 
solution  of  the  diffusion  equations  was  developed.  It  was  shown  that  the  process  model 
needs  further  improvement. 
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Figure  6.1:  Schematic  of  the  experimental  setup. 
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Figure  6.2:  Ten  replicate  experiments  conducted  on  LaNi5  electrode  at  open  circuit 
potential  under  potentiostatic  regulation. 
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Figure  6.3:  Plot  of  the  open  circuit  potential  as  a function  of  run  number  before  and  after 
each  scan  shown  in  Figure  6.2. 
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Figure  6.4:  Ten  replicate  experiments  conducted  on  LaNi5  electrode  zero  current  under 
galvanostatic  regulation. 
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Figure  6.5:  Plot  of  the  open  circuit  potential  as  a function  of  run  number  before  and  after 
each  scan  shown  in  Figure  6.4. 
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Figure  6.6:  Ten  replicate  experiments  conducted  on  LaNi5  electrode  at  a cathodic 
potential  of  -0.8V  (Hg/HgO). 
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Figure  6.7:  Ten  replicate  experiments  conducted  on  MmNi5  electrode  at  open  circuit 
potential  under  potentiostatic  regulation. 
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Figure  6.8:  Plot  of  the  open  circuit  potential  as  a function  of  run  number  before  and  after 
each  scan  shown  in  Figure  6.7. 
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Figure  6.9:  Ten  replicate  experiments  conducted  on  MmNi5  electrode  at  a cathodic 
potential  of  -1.0 V (Hg/HgO). 
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Figure  6.10:  Ten  replicate  experiments  conducted  on  MmNi5  electrode  at  an  anodic 
potential  of  -0.5  V (Hg/HgO). 
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Figure  6.11:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.2. 
The  circles  represent  the  experimental  data  and  the  solid  line  represents  complex  fit  to  the 
data. 
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Figure  6. 12:  Normalized  residuals  for  the  regression  shown  in  Figure  6. 1 1. 
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Figure  6.13:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.2. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the  real 
part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the  prediction. 
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Figure  6.14:  Residual  errors  for  the  regression  presented  in  Figure  6.13.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.15:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.2. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the 
imaginary  part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the 
prediction. 
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Figure  6.16:  Residual  errors  for  the  regression  presented  in  Figure  6.15.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.17:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.4. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the  real 
part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the  prediction. 
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Figure  6.18:  Residual  errors  for  the  regression  presented  in  Figure  6.17.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.19:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.4. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the 
imaginary  part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the 
prediction. 
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Figure  6.20:  Residual  errors  for  the  regression  presented  in  Figure  6.19.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.21:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.6. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the  real 
part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the  prediction. 
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Figure  6.22:  Residual  errors  for  the  regression  presented  in  Figure  6.21.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.23:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.6. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the 
imaginary  part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the 
prediction. 


172 


Figure  6.24:  Residual  errors  for  the  regression  presented  in  Figure  6.23.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.25:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.2.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.26:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.4.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.27:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.6.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.28:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.2.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.29:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.4.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.30:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.6.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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(Real  Fit) 


Figure  6.31:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.7. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the  real 
part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the  prediction. 
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Figure  6.32:  Residual  errors  for  the  regression  presented  in  Figure  6.31.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.33:  Results  of  regression  of  a measurement  model  to  data  shown  in  Figure  6.7. 
The  circles  represent  the  experimental  data  and  the  middle  line  represents  fit  to  the 
imaginary  part  of  the  data.  The  outer  lines  represent  the  95.4%  confidence  interval  for  the 
prediction. 
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Figure  6.34:  Residual  errors  for  the  regression  presented  in  Figure  6.33.  The  solid  lines 
represent  the  95.4%  confidence  interval. 
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Figure  6.35:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.7.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.36:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.9.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.37:  Plot  of  the  unfiltered  standard  deviation  for  the  data  shown  in  Figure  6.10. 
The  circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.38:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.7.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.39:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.9.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.40:  Plot  of  the  filtered  standard  deviation  for  the  data  shown  in  Figure  6.10.  The 
circles  represent  the  real  part  of  the  standard  deviation  and  the  triangles  represent  the 
imaginary  part. 
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Figure  6.41:  Results  of  the  regression  of  a process  model  (solid  line)  to  impedance  data 
(o)  for  LaNis  at  open  circuit  potential  under  potentiostatic  regulation.  No  weighting  was 
used  in  the  regression. 
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Figure  6.43:  Results  of  the  regression  of  a process  model  (solid  line)  to  impedance  data 
(o)  for  LaNi5  at  open  circuit  potential  under  potentiostatic  regulation.  Error  structure 
weighting  was  used  in  the  regression. 


Figure  6.44:  Normalized  residual  errors  for  the  regression  shown  in  Figure  6.43. 


CHAPTER  7 


MEASUREMENT  MODELS:  APPLICATION  TO  ELECTRO- 
HYDRODYNAMICAL  IMPEDANCE  SPECTROSCOPY 

In  the  previous  chapters,  the  concept  of  the  measurement  model  was  developed  and 
applied  to  electrochemical  impedance  spectroscopy.  In  this  chapter  the  measurement 
model  concept  is  extended  to  electro-hydrodynamical  impedance  spectroscopy  (EHD).  In 
this  chapter,  it  is  shown  that  identification  of  the  error  structure  for  EHD  and  regression  of 
a process  model  with  error  structure  weighting  eliminates  ambiguity  in  the  determination 
of  the  values  of  the  physical  parameters. 

7 . 1 Principles  of  EHD 

In  the  usual  application  of  electrochemical  impedance  spectroscopy,  a complex 
impedance  is  calculated  as  the  ratio  of  potential  to  current  under  small  perturbation  of 
current  (galvanostatic  regulation)  or  potential  (potential  regulation).  The  impedance  is 
measured  as  a function  of  the  frequency  of  the  perturbation,  and  regression  of  models  to 
the  resulting  spectra  yields  values  for  physical  properties.  In  recent  years  generalized 
impedance  techniques  have  been  introduced  in  which  a non-electrical  quantity  such  as 
pressure,  temperature,  magnetic  field,  and  light  intensity,  is  modulated  to  give  a current  or 
potential  response  [77,  78],  EHD  is  one  such  generalized  impedance  technique  in  which 
sinusoidal  modulation  of  disk  rotation  rate  drives  a sinusoidal  current  or  potential.  The 
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technique  has  been  applied  to  surface  or  electrode  processes  that  are  under  mass  transport 
control,  and  it  has  used  to  measure  diffusion  coefficients  of  ionic  species.  The  ionic 
diffusivities  are  usually  measured  by  measuring  the  Schmidt  number  (Sc  = v/D). 

The  advantage  of  using  EHD  over  usual  EIS  for  measuring  the  Sc  number  (or 
ionic  diffusivities)  is  illustrated  in  Figure  7.1.  A typical  current-potential  curve  for  a 
ferri/ferrocyanide  system  at  various  rotation  speeds  is  shown  in  Figure  7.1  [79],  With  EIS 
the  electrode  potential  is  modulated  sinusoidally  at  constant  disc  speed  and  the  Schmidt 
number  is  obtained  from  the  frequency  dependence  of  the  impedance  [80,  81], 
Examination  of  Figure  7.1  shows  that  to  obtain  a detectable  current,  the  ac  measurement 
must  be  made  on  the  slope  below  the  limiting  current  plateau.  At  currents  below  the 
limiting  current  plateau,  the  current  distribution  is  non-uniform  and  the  concentration  and 
surface  overpotential  vary  over  the  surface  of  the  electrode.  Therefore,  the  ac  method  is 
not  a general  reliable  method  for  measuring  the  Schmidt  number. 

The  EHD  method  in  contrast  uses  a sinusoidal  perturbation  in  the  speed  of  the 
electrode  at  a constant  electrode  potential  at  mass-transfer  limited  current.  The  resulting 
sinusoidal  current  is  measured  and  the  ratio  of  the  complex  current  and  the  amplitude  of 
the  disc  rotation  gives  the  EHD  transfer  function. 

The  EHD  transfer  function  has  been  related  to  the  EIS  transfer  function  [8 1],  For  a 
rotating  disc  electrode  the  instantaneous  state  of  the  electrochemical  system  can  be 
characterized  by  an  implicit  functional  relationship  such  as: 


(7.1) 
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where  I(t),  E(t)  and  Q(t)  are  the  instantaneous  values  of  current,  potential  and  rotation 
speed  respectively.  Differentiating  equation  7. 1 for  each  frequency  and  considering  that  I, 
E,  and  D.  are  independent  variables  two  by  two,  one  obtains, 
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is  now  defined  as  the  electro-hydrodynamical  impedance  under  galvanostatic  control  and 
is  denoted  as  Zehd.g,  and 
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is  defined  as  the  electro-hydrodynamical  impedance  under  potentiostatic  control  and 
denoted  <is  Zehdjp* 
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7.2  EHD  Experiments 

A series  of  EHD  experiments  were  conducted  on  a rotating  platinum  disc 
electrode.  The  experiments  were  designed  to  obtain  experimental  data  at  a wide  range  of 
experimental  conditions  so  that  the  data  could  be  analyzed  by  the  measurement  models  to 
determine  the  error  structure  of  the  data. 

The  disc  was  rotated  by  a rotating  disc  apparatus  developed  at  the  CNRS  [19]. 
The  potentials  and  currents  were  measured  and  controlled  by  a potentiostat.  A SI1250 
FRA  was  used  to  apply  the  sinusoidal  perturbation  and  calculate  the  transfer  function.  In 
addition  a low  pass  digital  filter  was  used  to  filter  the  input  and  output  to  and  from  the 
FRA.  Experiments  were  performed  with  the  filter  being  on-line  and  off-line  to  explore  its 
influence  on  the  noise. 

A ferri/ferro  cyanide  redox  couple  in  KC1  supporting  electrolyte  was  used  as  the 
experimental  system.  The  redox  reaction  for  this  system  is  very  fast  compared  to  the  mass 
transfer,  therefore,  this  system  is  suitable  for  the  comparison  of  diffusion  measurements. 
Equimolar  concentrations  of  potassium  ferricyanide  and  potassium  ferrocyanide  were 
dissolved  in  1M  potassium  chloride  solution.  The  different  concentrations  of  the  reacting 
species  used  were  0.1M,  0.0 1M,  and  0.00 1M.  Two  types  of  surface  treatments  were  used 
to  condition  the  working  electrode  before  the  start  of  the  experiment.  A set  of  experiments 
were  performed  on  an  electrode  which  was  pretreated  by  heating  it  to  a temperature  of 
100°C  for  two  hours.  Another  set  of  experiments  were  performed  where  the  electrode  was 
pretreated  by  cycling  the  potential  between  the  anodic  and  cathodic  potentials.  The 
experiments  were  conducted  at  two  mean  rpm’s  of  300  and  1200.  The  frequency  range  of 
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0.1Hz  to  100Hz  was  used  and  the  experiments  were  performed  at  a polarization  potential 
of -0.4V  (vs.  SCE).  A software  developed  in-house  at  CNRS  by  H.  Takenouti  was  used  to 
control  the  experiments  and  collect  the  data.  Twelve  replicate  runs  were  performed  at 
each  experimental  condition.  The  various  experimental  conditions  used  are  summarized  in 
Table  7.1.  A typical  set  of  EHD  data  is  shown  in  Figure  7.2.  The  experimental  conditions 
used  to  collect  the  data  shown  in  Figure  7.2  are  tabulated  in  Table  7.2. 

7.3  Measurement  Model  Analysis 

In  this  work,  the  measurement  model  concept  was  shown  to  apply  to  electro- 
hydrodynamic impedance  spectroscopy.  The  measurement  model  based  on  the  Voigt 
circuit  elements  shown  in  Figure  3.1  and  described  by  equation  3.1  was  regressed  to  the 
EHD  data  shown  in  Figure  7.2.  The  results  of  the  regression  are  shown  in  Figure  7.3  and 
the  examination  of  the  results  show  that  the  measurement  model  is  applicable  to  the  EHD 
data  with  the  term  R«  in  equation  3.1  set  to  zero. 

The  stochastic  components  of  the  errors  for  the  EHD  data  were  determined  using 
the  procedure  described  in  chapter  4.  The  results  of  the  analysis  for  the  data  shown  in 
Figure  7.4.  The  squares  in  Figure  7.4  represent  the  absolute  value  of  the  electro- 
hydrodynamical  impedance  for  the  data  shown  in  Figure  7.2.  The  circle  and  triangles 
represent  the  real  part  and  the  imaginary  part  of  the  filtered  standard  deviation.  In  contrast 
to  the  results  obtained  for  EIS  under  potentiostatic  control,  the  stochastic  noise  of  the 
EHD  measurements  appear , to  a first  approximation,  to  be  independent  of  frequency.  The 
standard  deviation  of  the  real  and  imaginary  parts  of  the  noise  are  equal,  a result  that  is 
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consistent  with  observations  made  for  EIS.  The  preliminary  model  for  the  error  structure 
for  EHD  data  was  determined  to  be 

ar  = =P\Zr\+8  (7.8) 

where  ov,  qj,  are  the  real  and  imaginary  part  of  the  standard  deviation,  Zr  is  the  real  part  of 
the  EHD  transfer  function  and  P and  5 are  parameters  which  depend  on  the 
instrumentation  being  used.  The  results  presented  in  Figure  7.4  also  illustrate  one  of  the 
problems  observed  during  the  regression  of  models  to  EHD  data,  i.e.,  that  the  noise  level 
becomes  a significant  fraction  of  the  signal  at  higher  frequencies. 

7.4  Process  Model  for  EHD 

The  electro-hydrodynamical  impedance  experiments  are  performed  by  modulating 
the  angular  velocity  Q of  the  rotating  disc  electrode  with  a fixed  amplitude  of  perturbation 
AQ  around  a mean  angular  velocity  Q0.  The  perturbation  frequency  © is  varied  to 
obtained  the  widest  possible  range.  The  angular  velocity  can  thus  be  written  as 

H = Q0  + ADcos(zyf)  (7  9) 

The  resulting  current  i,  at  fixed  potential,  can  be  expressed  as 

i = z'o  + A i cos  (co  t + <f>)  (7  10) 

where  z0is  the  average  current  and  M is  the  amplitude  of  the  sinusoidal  current  response 
and  <}>  is  the  phase  shift.  These  equations  can  then  be  introduced  into  the  differential 
equations  that  describe  the  physics  and  chemistry  of  the  system  to  calculate  the  transfer 
function  for  the  system  explicitly.  Tribollet  and  Newman  [82]  have  developed  the 
theoretical  transfer  function  for  the  system  which  is  of  the  form 
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where  /is  a function  they  tabulate,  A0  is  a constant  coefficient,  and  Sc  is  the  Schmidt 
number.  The  function^/ includes  a correction  for  the  diffusion  layer  not  being  infinitesimal, 
/.<?.,  for  the  Schmidt  number  not  being  infinite.  Equation  (7)  is  regressed  to  the 
experimental  data  to  obtain  the  values  for,40  and  Sc. 


7.5  Regression  of  Process  Model  to  Data 
Due  to  lack  of  information  on  the  error  structure,  modulus  weighting  has,  in  the 
past,  been  used  to  regress  the  EHD  model  to  the  data.  The  use  of  modulus  weighting 
implies  that  the  stochastic  noise  is  proportional  to  the  modulus  of  the  impedance,  an 
assumption  that  is  shown  in  Figure  7.4  to  be  incorrect. 

Inclusion  of  noisy  high  frequency  data  in  regression  under  modulus  weighting  has 
led  to  unrealistic  values  for  the  Schmidt  number.  The  problem  of  high  frequency  noise 
has,  in  the  past,  been  addressed  by  deleting  high  frequency  points.  Unfortunately, 
elimination  of  high  frequency  data  from  the  regression  leads  to  ambiguity  in  determination 
of  the  Schmidt  number.  The  results  of  the  regression  of  the  process  model  to  data  shown 
in  Figure  7.2,  as  a function  of  number  of  data  points  used  in  the  regression,  are  shown  in 
Figure  7.5.  The  absolute  value  of  the  impedance  was  used  to  weight  the  data.  Examination 
of  Figure  7.5  shows  that  although  the  estimate  for  cio  as  a function  of  number  data  points  is 
fairly  constant,  the  value  for  Sc  depends  on  the  number  of  data  points  included  in  the 
regression.  The  95%  confidence  interval  for  the  parameter  estimate  presented  in  Figure 
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7.5  as  upper  and  lower  solid  lines  also  reflects  large  uncertainty  in  the  estimate  for  Sc  and 
large  uncertainty  in  the  estimate  of  ao  when  more  than  30  points  are  used  in  the  regression. 

Another  approach  has  been  to  estimate  ao  independently  and  fix  it  during  the 
regression,  letting  Sc  be  the  only  unknown  parameter.  The  results  of  the  regression  are 
shown  in  Figure  7.6  where  a0  was  fixed  at  0.033481.  Examination  of  Figure  7.6  shows  a 
large  uncertainty  in  the  estimation  of  the  Sc  number  with  poor  confidence  in  the  estimate. 
In  fact  if  the  estimate  for  ao  is  changed  by  about  5%  there  is  a large  change  in  the  estimate 
of  Sc.  The  results  of  the  regression  with  a0  fixed  at  0.035  are  shown  in  Figure  7.7.  Hence 
to  get  a good  estimate  for  Sc  it  is  crucial  to  get  a good  estimate  for  a0  also.  Results  of  the 
regression  of  the  model  to  the  data,  with  proportional  weighting,  shown  in  Figure  7.8 
exhibit  similar  problems.  Use  of  the  measured  error  structure  for  the  EHD  measurements 
in  the  regression  given  by  equation  7.8  eliminates  ambiguity  in  the  identification  of  the 
Schmidt  number.  The  results  are  shown  in  Figure  7.9.  The  examination  of  the  results  show 
that  the  ambiguity  in  the  estimation  of  Sc  has  been  eliminated  except  when  too  few  points 
are  used  in  the  regression  resulting  in  the  loss  of  information.  The  confidence  interval  for 
the  estimates  is  also  seen  to  improve.  Use  of  the  error  structure  weighting  yields  an 
unambiguous  value  for  the  Schmidt  number  of  1383±15  as  compared  to  values  ranging 
from  1230  to  1460  when  modulus  weighting  is  employed  (see  Figure  7.2).  It  should  be 
noted  that  the  regressed  values  for  the  zero-frequency  asymptote  of  the  EHD  impedance 
shows  similar  ambiguity  when  modulus  weighting  is  employed  and  that  an  unambiguous 
value  is  obtained  when  the  error  structure  is  used  in  the  regression.  The  results  of  the 
regression  of  the  model  to  the  data  when  the  filter  was  on-line  are  shown  in  Figures  7.10- 
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7.12.  The  regression  results  exhibit  similar  features  and  it  is  observed  that  weighting  the 
regression  with  the  model  for  the  error  structure  gives  unambiguous  estimates  for  Sc. 

7.6  Conclusions 

Electro-hydrodynamical  impedance  spectroscopy  is  a useful  technique  for  the 
characterization  of  electrochemical  systems  at  mass  transfer  limited  currents,  but 
development  and  regression  of  appropriate  process  models  requires  characterization  of  the 
measurements  themselves.  The  measurement  model  provides  a means  of  characterizing 
the  measurement  characteristics  for  electro-hydrodynamical  impedance  spectroscopy  as 
well  as  electrochemical  impedance  spectroscopy.  Regression  employing  knowledge  of  the 
error  structure  yields  unambiguous  values  for  physical  properties  such  as  the  Schmidt 
number. 


Table  7,1:  Summary  of  the  experimental  conditions  used  to  collect  EHD  data. 


Working  Electrode 

Platinum  Disc 

Supporting  Electrolyte 

1MKC1 

Reacting  Species 

[Fe(CN)6f 

[Fe(CN)6;f 

Concentration  of  Reacting 
Species 

0.1M,  0.01M, 
0.00 1M 

RPM 

300,  1200 

Low  Pass  Digital  Filter 

On,  Off 

Frequency  Range 

O.lHzto  100Hz 

Polarization  Potential 

-0.4V  (SCE) 
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Working  Electrode 

Platinum  Disc 

Supporting  Electrolyte 

1MKC1 

Reacting  Species 

[Fe(CN)6f 

[Fe(CN)614~ 

Concentration  of  Reacting 
Species 

0.001M 

RPM 

300 

Low  Pass  Digital  Filter 

Off 

Frequency  Range 

O.lHzto  100Hz 

Polarization  Potential 

-0.4V  (SCE) 
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Figure  7.1:  Measured  electrochemical  current  vs  overpotential  for  five  values  of  disc  speed 
[79]  (reproduced  with  permission  from  the  publisher). 


205 


Figure  7.2.  Nyquist  plot  for  12  replicate  electro-hydrodynamical  impedance  experiments. 
The  experimental  parameters  used  are  shown  in  Table  7.2. 
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ZEHD,r 


Figure  7.3:  Results  of  the  regression  of  a measurement  model  to  a data  set  shown  in 
Figure  7.2.  The  circles  represent  the  experimental  data  and  the  solid  line  represents  the 
model. 
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Frequency,  Hz 


Figure  7.4:  The  error  structure  identified  for  EHD  measurements.  The  squares  represent 
the  absolute  value  of  electro-hydrodynamical  impedance  for  12  replicate  runs.  The  circles 
and  triangles  represent  the  real  part  and  imaginary  part  of  the  filtered  standard  deviation. 
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Frequency,  Hz 

0.10  1.00  10.0  100.0 


Frequency,  Hz 

0.10  1.00  10.0  100.0 


Figure  7.5:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression  to 
data  shown  in  Figure  7.2.  The  absolute  value  of  impedance  was  used  to  weight  the 
regression.  The  lines  represent  the  95%  confidence  interval  for  the  regressed  parameter. 
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Frequency,  Hz 

0.10  1.00  10.0  100.0 


Figure  7.6:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression.  The 
absolute  value  of  impedance  was  used  to  weight  the  regression.  The  value  of  a0  was  fixed 
at  0.033481.  The  lines  represent  the  95%  confidence  interval  for  the  regressed  parameter. 
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o Frequency,  Hz 

0.10  1.00  10.0  100.0 


Figure  7.7:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression.  The 
absolute  value  of  impedance  was  used  to  weight  the  regression.  The  value  of  ao  was  fixed 
at  0.035.  The  lines  represent  the  95%  confidence  interval  for  the  regressed  parameter. 
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Frequency,  Hz 

0.10  1.00  10.0  100.0 


Frequency,  Hz 

0.10  1.00  10.0  100.0 


Figure  7.8:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression. 
Proportional  weighting  was  used  in  the  regression.  The  lines  represent  the  95%  confidence 
interval  for  the  regressed  parameter. 
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Figure  7.9:  Schmidt  as  a function  of  number  of  data  points  used  in  the  regression.  The 
model  developed  in  this  work  for  the  standard  deviation  of  the  stochastic  noise  was  used 
to  weight  the  regression.  The  lines  represent  the  95%  confidence  interval  for  the 
regressed  parameter. 


213 


Frequency,  Hz 
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Frequency,  Hz 

0.10  1.00  10.0  100.0 


Figure  7.10:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression  to 
data  with  the  digital  filter-on.  The  absolute  value  of  impedance  was  used  to  weight  the 
regression.  The  lines  represent  the  95%  confidence  interval  for  the  regressed  parameter. 
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Frequency,  Hz 
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Frequency,  Hz 
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Figure  7.11:  Schmidt  as  a function  of  the  number  of  data  points  used  in  the  regression  to 
data  with  the  digital  filter  on.  Proportional  weighting  was  used  in  the  regression.  The 
lines  represent  the  95%  confidence  interval  for  the  regressed  parameter. 
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Figure  7.12:  Schmidt  as  a function  of  number  of  data  points  used  in  the  regression  to  data 
with  the  digital  filter  on.  The  model  developed  in  this  work  for  the  standard  deviation  of 
the  stochastic  noise  was  used  to  weight  the  regression.  The  lines  represent  the  95% 
confidence  interval  for  the  regressed  parameter. 


CHAPTER  8 


CONCLUSIONS 

A new  framework  for  model  identification  for  impedance  spectroscopy  was  proposed 
in  this  thesis.  The  approach  emphasizes:  use  of  deterministic  models  based  on  solution  of  the 
equations  governing  the  physics  and  chemistry  of  a given  system,  application  of  experimental 
techniques  sufficient  to  confirm  hypothesized  mechanisms,  and  application  of  statistical  tools  to 
ensure  that  measurement  characteristics  can  be  taken  into  account  when  comparing  models  to 
experimental  data. 

The  framework  of  the  general  approach  to  impedance  spectroscopy  was  divided  into 
three  main  categories: 

1.  Optimal  design  of  experiments 

2. Characterization  of  measurements 

3.  Model  identification 

An  in-house  software  to  run  impedance  experiments  optimally  was  developed.  The 
concept  of  a measurement  model  was  developed  and  was  used  to  characterize  the  stochastic 
and  systematic  errors.  An  empirical  model  for  the  stochastic  errors  was  proposed.  An 
algorithm  to  check  for  consistency  of  impedance  data  was  proposed.  These  concepts  were 
applied  to  a wide  class  of  systems.  It  was  shown  that  impedance  spectroscopy  is  an  effective 
technique  for  characterizing  metal  hydride  electrodes. 
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The  approach  was  also  applied  to  electro-hydrodynamical  impedance  of  a rotating 
platinum  disc  electrode.  It  was  shown  that  the  approach  helped  improve  the  characterization 
of  the  system  by  improving  the  regression  of  the  process  model  to  experimental  data. 


CHAPTER  9 


SUGGESTIONS  FOR  FUTURE  WORK 

This  inferential  character  of  impedance  spectroscopy  and  the  ambiguity  inherent  in 
the  interpretation  of  impedance  spectra  give  rise  to  several  interesting  questions  during 
model  identification.  For  example: 

•Do  the  data  satisfy  the  Kramers-Kronig  transforms  [5,6]  (i.e.,  are  the  data 
representative  of  a stationary,  causal,  and  linear  system)? 

•What  is  the  “correct”  model  for  a data  set  (i.e.,  should  a data  set  be  interpreted  in 
terms  of  coupled  reactions,  surface  roughening,  frequency  dispersion,  or  something  else)? 

•Do  the  models  really  provide  a statistically  significant  fit  to  the  data? 

•What  portion  of  the  spectrum  can  be  interpreted  in  terms  of  a model  that  assumes 
stationary  behavior? 

•Is  the  experimental  protocol  optimized  with  respect  to  minimizing  errors  in  the  data? 

The  measurement  model  appears  to  provide  a useful  framework  for  a statistical 
analysis  of  impedance  data.  The  approach  is  currently  limited  to  evaluating  consistency  of 
data  with  the  Kramers-Kronig  relations,  to  filtering  data  for  lack  of  replicacy  (a  step  in  the 
experimental  identification  of  the  error  structure  of  impedance  measurements),  and  to 
providing  a rough  guide  for  model  identification.  Work  is  needed  to  improve  the  use  of 
measurement  models  for  model  identification,  for  characterization  of  industrial  data 
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(where  extrapolation  is  often  more  important  than  model  identification,  and  for 
characterization  of  non-stationary  phenomena.  The  model  for  the  error  structure  proposed 
in  this  work  is  preliminary  and  needs  to  be  improved. 

The  application  of  impedance  spectroscopy  to  metal  hydride  systems  requires  more 
experiments  at  a wide  range  of  experimental  conditions,  e.g.,  different  state  of  charge.  A 
more  detailed  process  model  needs  to  be  developed  to  characterize  the  complete  systems. 
Some  of  the  processes  which  need  to  be  incorporated  into  the  model  are  the  corrosion 
reaction,  potential  dependence  of  the  diffusion  coefficient,  lattice  expansion,  and  the 
formation  of  the  hydride  etc. 


APPENDIX  A 


FORTRAN  CODE  FOR  MEASUREMENT  MODEL 

This  program  estimates  the  parameters  for  a Voigt  Circuit  (see  equation  2.1)  from 
impedance  data.  The  program  is  based  on  the  original  code  developed  by  Prof.  Luis  H. 
Garcia-Rubio  for  the  estimation  of  parameter  estimates  for  Debye  oscillator  circuits  from 
optical  spectroscopy  data.  The  regression  algorithm  developed  and  implemented  by  Prof. 
Garcia-Rubio  is  not  included  here  but  is  used  in  its  original  form  to  perform  the 
regressions. 

The  program  allows  the  user  the  following  weighting  choices: 

1.  no  weighting 

2.  weighting  by  error  structure  proposed  in  this  work 

3.  proportional  weighting 

4.  modulus  weighting 

The  program  was  implemented  on  a Dec3 100  workstation  with  an  Ultrix  operating 
system.  The  program  can  handle  a maximum  of  10  data  files  with  a maximum  100  data 
points  in  each  file.  The  maximum  number  of  Voigt  elements  allowed  is  15. 
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c 

c 

c c 

c c 

c This  program  estimates  the  parameters  of  a Voigt  Circuit  c 
c from  impedance  and  resistance  data  c 

c c 

c c 

c c 

c Pankaj  Agarwal  and  Mark  E.  Orazem  November/1994 

c c 

c c 

c data  epsl,eps2/1.0e-08,1.0e-03/ 
c scrat=5*np+np**2+2*nob+np*nob 


c fitype  = C for  complex  I for  imaginary  R for  real  fit 
c nda  = no  of  data  points  in  each  dataset 
c nfile  = no  of  files  to  be  used  in  the  regression 
c nrand  = no  of  exp  for  Monte  Carlo 
c ntest  = the  no  of  the  frequency  for  which  data  is  to  stored 
c after  Monte  Carlo 

c nweight  = 0 no  weighting 

c nweight  = 1 error  structre  used  as  weights  ( data  values) 

c nweight  = -1  error  structre  used  as  weights  (function  values) 

c nweight  = 2 weighting  equal  to  data  or  mean  value  of  data  (for 

c multiple  data  sets) 

c nweight  = -2  weighting  equal  to  function  value 

c nweight  = 3 absolute  value  of  impedance  ( data) 

c nweight  = -3  absolute  value  of  impedance  (function) 

c nweight  = 4 actual  variance  used  as  weights 

c nweight  = 5 then  the  weights  are  read  from  the  file  data,  weight, 

c per=  the  percent  error  to  estimate  the  noise  level  in  the  data, 
c nose  = no  of  oscillators  max=14 
c epsl,eps2  = tolerances  for  the  regression 

c mit  = max  no  of  iterations 
c ndat  = total  no  of  data  points  (nfile*nda) 


c 


c 


c c 

c 


common/osc/nosc,tau,eo,eou,eol 

common/sig/nsign,nsignsol,nsol 

common/erst/alpha,beta,gamma,delta 

common/serr/nalpha,nbeta,ngamma,ndelta 

common/mrs/vecres 

common/opt/nda,nfile,nrand,ntest,nweight,mit 

common/perc/per 

common/type/fitype,  comment 

common/var/datarm,dataim,datars,datais 

common/cl/refr,refi,wav 

common/weight/weightr,weighti 

common/calcf7funr,funi 

common/par/bpar,upar 

common/modelc/zrm,zim 
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common/err/er,ei 

common/frequency/fre 

common/stat/ssqr,ssqi,ssq,erssq,eissq,essq,ssdatr,ssdati,ssdat 

common/eps/eps  l,eps2 

common/sum  sq/sumb 

common/erul/zru,zrl,ziu,zil 

common/ratio/ratior,ratioi 

common/akaike/fpe,wfpe,aic,waic 

dimension  ratior(  10  l),ratioi(  10  l),vecres(10 1) 
dimension  zru(101),zrl(101),ziu(101),zil(101) 
dimension  weighti(101),weightr(101),bpar(34),upar(34) 
dimension  er(101),ei(101),funr(101),funi(101) 
dimension  y(1600),th(34),diff(34),signs(34),wave(1000) 
dimension  scrat(20000),fre(101),taul(4,15),tauu(4,15) 
dimension  nsign(4,15),tau(4,15) 
dimension  refr(1000),refi(1000),wav(1000) 
dimension  datarm(101),dataim(101),datars(101),datais(101) 
dimension  zrm(1000),zim(1000) 
character  fitype*l,comment*75 
c 

external  model 
c 

data  flam,fhu/.0 1,10.0/ 
c 

nsti  = 5 
nsto  = 6 
pi=3. 14 159265 
c 

open(  1 1 ,file='out',  status-unknown') 
open(15,file='out.new', status-unknown') 
open(16,file-out.err',status- unknown') 
open(  17,file='measure.  error',  status-unknown') 
open(20,file='test.  out', status-unknown') 
c 

c read  program  parameters  from  file  'input' 


call  readinput 
ndat=nda*nfile 

call  mergefile(wave,refr,refi,nfile,nda,nweight,fre) 

c mergefile  reads  the  data  and  calculates  the  mean  and  variance  for  the  c 

do  20  i=l,ndat 

refr(i)=refr(i) 
refi(i)=-l  .0*refi(i) 
wav(i)=wave(i)*2*pi 

npoint=mod(i,nda) 
if  (npoint.eq.O)  npoint=nda 


data 
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c 

c statiscally  weight  the  data  but  retain  sign 
c 

if  (fitype.eq."CM)  then 

y(i)  = refi(i)/abs(weighti(npoint)) 
y(i+ndat)  = refr(i)/abs(weightr(npoint)) 

endif 

if  (fitype.eq."I")  then 

y(i)  = refi(i)/abs(weighti(npoint)) 

endif 

if  (fitype.eq."R")  then 

y(i)  = refr(i)/abs(weightr(npoint)) 

endif 

c 

20  continue 

c 

c 

c 

c set  uwhaus  parameters  for  estimation 
c nob  = total  no  observations  eq  to  ndat  for  real  and  im  fit 
c eq  to  2*  ndat  for  complex  fit 


if  (fitype.eq.T)  eoi=eo 


c 

nob=ndat 

c 

if  (fitype.eq."C")  nob=ndat*2 
c 
c 

ii=l 

np=0 

do  35  i=l,nosc 

do  35  j=l,4 

if  (nsign(j,i).ne.O)  then 

if  (nsign(j,i).gt.O)  then 
signs(ii)=1.0 

else 

signs(ii)=-1.0 

endif 
th(ii)=tau(j,i) 
ii=ii+l 
np=np+l 
endif 
35  continue 
c 

if  (nsignsol.ne.O)  then 
np=np+l 
th(np)  = eo 


signs(np)=1.0 

endif 

nsol=np 

ierr=0 

if  (nalpha.ne.O)  then 
np=np+l 
ierr=ierr+l 
th(nsol+ierr)=alpha 
signs(nsol+ierr)=1.0 

if(nalpha.lt.O)  signs(nsol+ierr)=-1.0 

endif 

if  (nbeta.ne.O)  then 

np=np+l 
ierr=ierr+l 
th(nsol+ierr)=beta 
signs(nsol+ierr)=l  .0 

if(nbeta.lt.O)  signs(nsol+ierr)=-1.0 

endif 

if  (ngamma.ne.O)  then 
np=np+l 
ierr=ierr+l 
th(nsol+ierr)=gamma 
signs(nsol+ierr)=  1 . 0 

if(ngamma.lt.O)  signs(nsol+ierr)=-1.0 

endif 

if  (ndelta.ne.O)  then 

np=np+l 

ierr=ierr+l 

th(nsol+ierr)=delta 

signs(nsol+ierr)=1.0 

if(ngamxna.lt.O)  signs(nsol+ierr)=-1.0 

endif 

do  40  i=l,np 
diff(i)=0.01 
continue 


call  uwhaus  for  parameter  estimation 

call  uwhaus(nprob,model,nob,y,np,th,diff,signs,eps  l,eps2, 
1 mit,flam,fnu,scrat) 


retrieve  the  parameter  values 


if  (fitype.eq."C")  then 

print*,'  calculated  parameter  values  for  COMPLEX  fit ' 
endif 

if  (fitype.eq."R")  then 

print*,'calculated  parameter  values  for  REAL  fit ' 
endif 

if  (fitype.eq.'T')  then 

print*, 'calculated  parameter  values  for  IMAGINARY  fit ' 
endif 

ii=l 

do  36  i=l,nosc 
do  36  j=l,4 

if  (nsign(j,i).ne.O)  then 
tau(j,i)  = th(ii) 
tauu(j,i)=  upar(ii) 
taulO,i)=  bpar(ii) 
ii=ii+l 

endif 

continue 

if  (nsignsol.ne.O)  then 
eo=th(nsol) 
eou=upar(nsol) 
eol=bpar(nsol) 
endif 


ierr=0 

if  (nalpha.ne.O)  then 

ierr=ierr+l 

alpha=th(nsol+ierr) 

endif 

if  (nbeta.ne.O)  then 

ierr=ierr+l 

beta=th(nsol+ierr) 

endif 

if  (ngamma.ne.O)  then 
ierr=ierr+l 
gamma=th(nsol+ierr) 

endif 

if  (ndelta.ne.O)  then 

ierr=ierr+l 

delta=th(nsol+ierr) 


endif 


calculate  the  impedance  values 


ssqr=0.0 
ssqi=0.0 
do  50  i=l,nda 
do  51  k=l,nfile 
ncr=(k-l)*nda 


w = wav(i+ncr) 
call  debye(nosc,w,tau,xip,xipp) 
ep  = eo  + xip 
epp=  xipp 


if  (ncr.eq.O)  then 
zrm(i)=ep 
zim(i)=epp 
endif 

write(l  1,*)  wav(i+ncr),refr(i+ncr),ep,refi(i+ncr),epp 


this  section  has  added  on  12/13/92 

phid=atan(refi(i+ncr)/refr(i+ncr)) 

phic=atan(epp/ep) 

zabsd=  (reff(i+ncr)**2  + refi(i+ncr)**2)**0.5 

zabsc=  (ep**2+epp**2)**0.5 

errphi  = phid-phic 

errzabs  = zabsd-zabsc 

errzr  = refr(i+ncr)  -ep 

errzi  = refi(i+ncr)  -epp 

write(15,*)  wav(i+ncr),zabsd,zabsc,phid,phic 
write(16,*)  wav(i+ncr), errzr, errzi, errzabs, errphi 

new  section  ends 


ssqr  = ssqr  + (refr(i+ncr)-ep)**2 
ssqi  = ssqi  + (refi(i+ncr)-epp)**2 


continue 

continue 

erssq=0.0 

eissq=0.0 

ssdatr=0.0 

ssdati=0.0 
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do  52  i=l,nda 

erssq=erssq+er(i)**2 

eissq=eissq+ei(i)**2 

ssdatr=ssdatr+datarm(i)**2 

ssdati=ssdati+dataim(i)**2 


52  continue 


erssq=eissq/nda 

eissq=eissq/nda 

ssdatr=ssdatr/nda 

ssdati=ssdati/nda 


ssqr=ssqr/ndat 

ssqi=ssqi7ndat 

ssq  = ssqr  + ssqi 
essq  = erssq  +eissq 
ssdat=  ssdatr  + ssdati 
c 
c 


c calculte  the  Akaike  information  index 
if  (fitype.eq."C")  then 

fpe=ssq*(1.0+np/(nob*  1 ,0))/(1.0-np/(nob*  1 .0)) 
aic=log((1.0  + (2.0*np)/nob)*ssq) 
print*, 'no  of  data  pts  =',nob,'  no  of  parameters  =',np 
print*, 'fpe=',fpe,'  aic=',aic,'  sum  of  sq  - ,ssq 

endif 


if  (fitype.eq."R")  then 

fpe=ssqr*(1.0+np/(nob*  1 ,0))/(1.0-np/(nob*  1 .0)) 
aic=log((1.0  + (2.0*np)/nob)*ssqr) 
print*, 'no  of  data  pts  - ,nob,'  no  of  parameters  =',np 
print*, 'fpe=',fpe,'  aic=',aic,'  sum  of  sq  - ,ssqr 

endif 


if  (fitype.eq.'T')  then 

fpe=ssqi*(l  .0+np/(nob*  1 .0))/(1.0-np/(nob*  1.0)) 
aic=log((1.0  + (2.0*np)/nob)*ssqi) 
print*, 'no  of  data  pts  =',nob,'  no  of  parameters  =',np 
print*, 'fpe=’,fpe,'  aic=',aic,'  sum  of  sq  =',ssqi 


endif 


wfpe=sumb*(1.0+np/(nob*  1.0))/(1.0-np/(nob*  1.0)) 

waic=log((1.0  + (2.0*  np)/nob)  * sumb) 

print*, 'wfpe-,wfpe,'  waic=',waic,'  wssq=',sumb 


write(22,'(a)')  fitype 
write(22,*)  ndat 
write(22,*)  nose 


do  69  jj=l,np 

69  write(22,*)  upar(jj),bpar(ij) 


do  91  i=l,nosc 
print*,i 

print*,'  Rw=  ',tau(l,i),'  Tw=  ',tau(2,i) 
print*,'  R=  ’,tau(3,i),'  T=  \tau(4,i) 


91  continue 

print*,'  eo=  ',eo 

print*,'  real  sum  of  squares  res  err  = ',ssqr 
print*,'  imaginary  sum  of  squares  res  err  = ',ssqi 
print*,'  total  sum  of  squares  res  err  = ',ssq 

print*,'  real  sum  of  squares  err  = ',erssq 
print*,'  imaginary  sum  of  squares  err  = ',eissq 
print*,'  total  sum  of  squares  err  = ',essq 

print*,'  real  sum  of  squares  = ',ssdatr 
print*,'  imaginary  sum  of  squares  = ',ssdati 
print*,'  total  sum  of  squares  = ',ssdat 


call  writeoutput 
call  parfile(np) 
write(14,*)  ssqr,ssqi,ssq 
write(14,*)  erssq,eissq,essq 
write(14,*)  ssdatr,ssdati,ssdat 

open(3  l,file='ratio.file',status='unknown') 

do  199  i=l,nda 

write(17,*)  fre(i),zru(i),zrl(i),ziu(i),zil(i) 
write(31,*)  fre(i),ratior(i),ratioi(i) 

199  continue 


close(ll) 

close(15) 

close(16) 


close(17) 
close(20) 
close(21) 
close(22) 
close(3 1) 


print*, 'type  1 if  you  want  calculate  the  confidence  interval' 
read*,nconf 
if  (nconf.ne.  1)  goto  119 
c print*,fitype,nconf 

call  confidence(fitype,ndat,np-ierr) 
c print*,  fitype,nconf+l 

c 

119  continue 
stop 
end 

c 

c 

c 

c SUBROUTINE  MODEB03. FOR 

c 

c 

SUBROUTINE  MODEL(nprob,th,f,nob,np) 
c 

common/osc/nosc,tau,eo,eou,eol 

common/sig/nsign,nsignsol,nsol 

common/erst/alpha,beta,gamma,delta 

common/serr/nalpha,nbeta,ngamma,ndelta 

common/mrs/vecres 

common/opt/nda,nfile,nrand,ntest,nweight,mit 
common/type/fitype, comment 
common/var/datarm,dataim,datars,datais 
common/cl/reff,refi,wav 
common/weight/weightr,weighti 
common/calcfifiinr,funi 
c 

character  fitype*l,comment*75 
dimension  th(34),f(1600),nsign(4,15),tau(4,15) 
dimension  weighti(101),weightr(101),funr(101),funi(101) 
dimension  datarm(101),dataim(101),datars(101),datais(101) 
dimension  refr( 1000), refi( 1000),  wav( 1000),vecres(  10 1) 
c 

c retrieve  the  parameter  values 

c 


ii=l 

do  35  i=l,nosc 
do  35  j=l,4 

if  (nsign(j,i).ne.O)  then 
tau(j,i)=  th(ii) 


endif 


ii=ii+l 


35  continue 
c 

if  (nsignsol.ne.O)  then 
eo=th(nsol) 
endif 

ierr=0 

if(nalpha.ne.O)  then 

ierr=ierr+l 

alpha=th(nsol+ierr) 

endif 

if(nbeta.ne.O)  then 

ierr=ierr+l 

beta=th(nsol+ierr) 

endif 

if(ngamma.ne.O)  then 

ierr=ierr+l 

gamma=th(nsol+ierr) 

endif 

if(ndelta.ne.O)  then 

ierr=ierr+l 

delta=th(nsol+ierr) 

endif 


c calculate  the  impedance  and  the  resistivity 
c 

ndat  = nob 

if  (fitype.eq."C")  then 
ndat  = nob/2 

endif 

if  (nweight.lt.  0)  then 

do  100  i=l,ndat 

w = wav(i) 
npoint=mod(i,nda) 
if  (npoint.eq.0)  npoint=nda 
call  debye(nosc,w,tau,xip,xipp) 
ep  = eo  + xip 
epp=  xipp 
funr(npoint)=ep 
funi(npoint)=epp 


100 


continue 
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endif 

if  (abs(nweight).eq.  1)  call  measres(nweight,np,nda) 
if  (abs(nweight).eq.2)  call  statwt(nweight,np,nda) 
if  (abs(nweight).eq.3)  call  propwt(nweight,np,nda) 


do  153  i=l,ndat 

w = wav(i) 
npoint=mod(i,nda) 
if  (npoint.eq.O)  npoint=nda 
call  debye(nosc,w,tau,xip,xipp) 

ep  = eo  + xip 

epp=  xipp 

if  (fitype.eq."C")  then 

f(i)  = epp/abs(weighti(npoint)) 
f(i+ndat)  = ep/abs(weightr(npoint)) 

endif 

c 

c 

if  (fitype.eq.T')  then 

f(i)  = epp/abs(weighti(npoint)) 

endif 

c 

c 

if  (fitype.eq."R")  then 

f(i)  = ep/abs(weightr(npoint)) 

endif 

c 

153  continue 
return 
end 
c 

c subroutine  read  input  reads  the  input  parameters  from  file  input 
c 

c fitype  = C for  complex  I for  imaginary  R for  real  fit 
c nda  = no  of  data  points  in  each  dataset 
c nfile  = no  of  files  to  be  used  in  the  regression 
c nrand  = no  of  exp  for  Monte  Carlo 
c ntest  = the  no  of  the  frequency  for  which  data  is  to  stored 
c after  Monte  Carlo 

c nweight  = 0 no  weighting 

c nweight  = 1 error  structre  used  as  weights  ( data  values) 

c nweight  = -1  error  structre  used  as  weights  (function  values) 

c nweight  = 2 weighting  equal  to  data  or  mean  value  of  data  (for 

c multiple  data  sets) 

c nweight  = -2  weighting  equal  to  function  value 

c nweight  = 3 absolute  value  of  impedance  ( data) 

c nweight  = -3  absolute  value  of  impedance  (function) 

c nweight  = 4 actual  variance  used  as  weights 

c nweight  = 5 then  the  weights  are  read  from  the  file  data. weight. 
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c per=  the  percent  error  to  estimate  the  noise  level  in  the  data, 
c nose  = no  of  oscillators  max=14 
c epsl,eps2  = tolerances  for  the  regression 
c mit  = max  no  of  iterations 
c ndat  = total  no  of  data  points  (nfile*nda) 
c 


SUBROUTINE  readinput 

common/osc/nosc,tau,eo,eou,eol 

common/sig/nsign,nsignsol,nsol 

common/erst/alpha,beta,gamma,delta 

common/serr/nalpha,nbeta,ngamma,ndelta 

common/perc/per 

common/opt/nda,nfile,nrand,ntest,nweight,mit 
common/type/fitype, comment 
common/eps/epsl,eps2 

character  fitype*  l,comment*75 
dimension  tau(4, 15),nsign(4, 15) 

open(13,file=,input,,status='old') 

read(13,'(a)')  comment 
read(13,'(a)')  fitype 

read(13,*)  nda,nfile,nrand,ntest,nweight 

read(13,*)  per 

read(13,*)  nose 

read(13,*)  epsl,eps2 

read(13,*)  mit 

do  31  i=l,nosc 

read(  13 , *)  tau(  l,i),tau(2,i),tau(3  ,i),tau(4,i) 
read(  13 , *)  nsign(  1 ,i),nsign(2,i),nsign(3  ,i),nsign(4,i) 
3 1 continue 

if  (fitype.  ne.  "I")  then 
read(13,*)  eo 
read(13,*)  nsignsol 
endif 

if  (fitype.  eq.  "I")  then 
read(13,*)  eo,eou,eol 
read(13,*)  nsignsol 
endif 

read(13,*)  alpha,beta,gamma,delta 
read(13,*)  nalpha,nbeta,ngamma,ndelta 

close(13) 

return 

end 


SUBROUTINE  writeoutput 


common/osc/nosc,tau,eo,eou,eol 
common/  sig/nsign,  nsignsol,  nsol 
common/erst/alpha,beta,gamma,delta 
common/serr/nalpha,nbeta,ngamma,ndelta 
common/perc/per 

common/opt/nda,nfile,nrand,ntest,nweight,mit 
common/type/fitype, comment 
common/eps/eps  1 ,eps2 

character  fitype*  l,comment*75 
dimension  tau(4,15),nsign(4,15) 

open(14,file='output',status='unknown') 

write(14,'(a)')  comment 
write(14,'(a)')  fitype 

write/14,*)  nda,nfile,nrand,ntest,nweight 
write/ 14,*)  per 
write/ 14,*)  nose 
write(14,*)  epsl,eps2 
write(14,*)  mit 

do  31  i=l,nosc 

write/ 14,  *)  tau(  1 ,i),tau(2,i),tau(3  ,i),tau(4,i) 
write/ 14, *)  nsign/ 1 ,i),nsign(2,i),nsign(3 ,i),nsign(4,i) 
continue 

write/14,*)  eo,eou,eol 
write/14,*)  nsignsol 
write/14,*)  alpha,beta,gamma,delta 
write(14,*)  nalpha,nbeta,ngamma,ndelta 


close(14) 

return 

end 

subroutine  mergefile(wave,refr,refi,nfile,nda,nweight,fie) 


common/var/datarm,dataim,datars,datais 

common/weight/weightr,weighti 

common/files/farr 

common/calcf7funr,funi 

common/mrs/vecres 

dimension  refr/1000),refi(1000),wave(1000),fre(101) 
dimension  datarm(101),dataim(101),datars(101),datais(101) 
dimension  datar(10,101),datai(10,101),vecres(101) 


dimension  weighti(10  l),weightr(10  l),funr(101),funi(  10 1) 
character  fiiame*40,farr(10)*40 
logical  fexist*4 
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nsti  = 5 
nsto  = 6 


c read  the  data  files 
do  11  i=l,nfile 

10  write(nsto,*)  'enter  filename  for  impedance  data' 
read(nsti,'(A)')  fname 
farr(i)  = fname 

inquire(file=fhame,exist=fexist) 
if(.not.fexist)  then 

write(nsto,*)  fname, 'does  not  exist,enter  new  name' 
go  to  10 

endif 

open(25,file=fhame, status- old') 
do  21  k=l,nda 

nc=(i-l)*nda 

read(25,*)  wave(nc+k),refr(nc+k),refi(nc+k), 

1 garb,garb,res 

datar(i,k)=reff(nc+k) 
datai(i,k)=refi(nc+k) 
if  (nc.eq.O)  then 

fre(k)=wave(k) 

if  (res.gt.  1.)  vecres(k)=10**(((res-2)*10)-l) 
if  (res.lt.  1.)  vecres(k)=10**((res*10)-l) 
endif 


21  continue 

close(25) 

11  continue 

c calculate  the  mean  and  the  standard  deviation  for  the  data 

do  31  ii=l,nda 

sumrm=0.0 

sumim=0.0 

do  41  kk=l,nfile 

sumrm=sumnn+datar(kk,ii) 

sumim=sumim+datai(kk,ii) 

41  continue 
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datarm(ii)=sumrin/nfile 

dataim(ii)=sumim/nfile 


sumrs=0.0 

sumis=0.0 

do  51  jj=l,nfile 

sumrs=sumrs+(datar(jj,ii)-datarm(ii))**2 

sumis=sumis+(datai(ij,ii)-dataim(ii))**2 

51  continue 

datars(ii)=(sumrs/nfile)**0.5 
datais(ii)=(sumis/nfile)*  *0. 5 

3 1 continue 


if  (nweight.eq.0)  then 
do  129  i=l,nda 

weightr(i)=1.0 

weighti(i)=1.0 

129  continue 
endif 


if  (nweight.lt.O)  then 

do  171  nc=l,nda 

funr(nc)=datarm(nc) 

funi(nc)=dataim(nc) 

171  continue 
endif 

if  (abs(nweight).eq.l)  call  measres(nweight,np,nda) 
if  (abs(nweight).eq.2)  call  statwt(nweight,np,nda) 
if  (abs(nweight).eq.3)  call  propwt(nweight,np,nda) 


if  (nweight.eq.4)  then 
do  149  i=l,nda 

weightr(i)=datars(i) 

weighti(i)=datais(i) 

149  continue 

endif 

if  (nweight.eq.5)  then 

open(30,file='data.weight',status='old') 
do  91  nv=l,nda 

read(30,*)  ffe(nv),weightr(nv),weighti(nv) 

91  continue 

close(30) 
endif 
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if  (nfile.eq.  1)  then 
do  92  jj=l,nda 

datarm(jj)=refr(ij) 

dataim(ij)=refi(ij) 

datars(ij)=refr(ij) 

datais(jj)=refi(ij) 

92  continue 
endif 

open(l,file='data.var', status- unknown') 
do  111  1=1, nda 

write(l,*)fre(l),datarm(l),datars(l),dataim(l),datais(l) 
111  continue 


close(l) 

return 

end 


subroutine  measres(nweight,np,nda) 

common/var/datarm,dataim,datars,datais 

common/err/er,ei 

common/weight/weightr,weighti 

common/mrs/vecres 

common/erst/alpha,beta,gamma,delta 

common/calcf/funr,funi 

common/frequency/fre 

common/erul/zru,zrl,ziu,zil 

conunon/ratio/ratior,ratioi 

dimension  fre(101),vecres(101) 
dimension  weighti(101),weightr(101) 
dimension  er(101),ei(101),funr(101),funi(101) 
dimension  datarm(101),dataim(101),datars(101),datais(101) 
dimension  zru(101),zrl(101),ziu(101),zil(101) 
dimension  em(101),ein(101),ratior(101),ratioi(101) 

if  (nweight.eq.-l)  then 

do  11  i=l,nda 

az  = ((funr(i)**2+funi(i)**2)**0.5) 

er(i)=  alpha*abs(funi(i))+beta*abs(funr(i))+ 

1 (gamma’"(az’,‘’''2))Arecres(i)  + delta 

ei(i)=  alpha*abs(funi(i))+beta*abs(funr(i))+ 

1 (gamma*  (az*  *2))/vecres(i)  + delta 

azn  = ((datarm(i)**2+dataim(i)**2)**0.5) 


em(i)=  alpha*abs(dataim(i))+beta*abs(datarm(i))+ 
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1 (gamma*  (azn*  *2))/vecres(i)  + delta 

ein(i)=  alpha*abs(dataim(i))+beta*abs(datarm(i))+ 
1 (gamma*(azn**2))/vecres(i)  + delta 

ratior(i)=  er(i)/em(i) 
ratioi(i)=  ei(i)/ein(i) 


11  continue 

do  139  i=l,nda 

weightr(i)=er(i) 

weighti(i)=ei(i) 

139  continue 
endif 

if  (nweight.eq.l)  then 

do  12  i=l,nda 

az  = ((datarm(i)**2+dataim(i)**2)**0.5) 

er(i)=  alpha*abs(dataim(i))+beta*abs(datarm(i))+ 
1 (gamma*(az**2))/vecres(i)  + delta 

ei(i)=  alpha*abs(dataim(i))+beta*abs(datarm(i))+ 
1 (gamma*(az**2))/vecres(i)  + delta 


12  continue 

do  140  i=l,nda 

weightr(i)=er(i) 

weighti(i)=ei(i) 

140  continue 
endif 


do  13  i=l,nda 


zru(i)=datarm(i)+er(i) 

zrl(i)=datarm(i)-er(i) 

ziu(i)=dataim(i)+ei(i) 

zil(i)=dataim(i)-ei(i) 


13  continue 


return 

end 
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subroutine  statwt(nweight,np,nda) 
common/var/datam,dataim,datars,datais 
common/weight/weightr,weighti 
common/calcf7funr,funi 

dimension  weighti(  10  l),weightr(  10  l),funr(10  l),funi(  101) 
dimension  datarm(101),dataim(101),datars(101),datais(101) 

if  (nweight.eq.-2)  then 
do  119  i=l,nda 

weightr(i)=funr(i) 

weighti(i)=funi(i) 

119  continue 
endif 

if  (nweight.eq.2)  then 
do  120  i=l,nda 

weightr(i)=datarm(i) 

weighti(i)=dataim(i) 

120  continue 

endif 

return 

end 


subroutine  propwt(nweight,np,nda) 
common/var/datarm,dataim,datars,datais 
common/weight/weightr,weighti 
common/cald7funr,funi 

dimension  weighti(101),weightr(101) 
dimension  fimr(101),funi(101) 

dimension  datarm(101),dataim(101),datars(101),datais(101) 
if  (nweight.eq.3)  then 

do  11  i=l,nda 

weightr(i)=(datarm(i)**2+dataim(i)**2)**0.5 

weighti(i)=(datarm(i)**2+dataim(i)**2)**0.5 

11  continue 

endif 


if  (nweight.eq.-3)  then 

do  12  i=l,nda 

weightr(i)=(funr(i)**2+funi(i)**2)**0.5 
weighti(i)=(funr(i)*  ^ *2+funi(i),,'  *2)*  *0.5 

12  continue 


nnn  non 
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endif 

end 

subroutine  deb(rd,td,zrd,zid,w) 

complex  cnum,cden,csum 

cnum  = cmplx(rd,0.0) 
cden  = cmplx(1.0,(-td*w)) 

csum  = cnum/cden 

zrd=real(csum) 

zid=aimag(csum) 

return 

end 


GENERALIZED  FINITE  WARBURG 

SUBROUTINE  CAPWarburg  (Frequency,  Rww,  Tww,Cww,Rwwc,Zrww,  Ziww) 
IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

COMPLEX*  16  CTEMP,EXPSAV, CTANH, DISTEL 
real*4  rd,td,zrd,zid,fre 


if  (rww.le.l.0d-10.or.tww.le.l.0d-10)  then 

rd=sngl(rwwc) 

td=sngl(cww) 

fre=sngl(frequency) 

call  deb(rd,td,zrd,zid,fre) 

zrww=dble(zrd) 

ziww=dble(zid) 

return 

endif 

phi=0.5 

PI=  3.141592654 
PI02  = 3.141592654/2.0 
ARG  = PI02*PHI 

CTEMP  = (Tww*frequency)**PHI  * DCMPLX(DCOS(ARG),  DSIN(ARG)) 

CALCULATE  COMPLEX  HYPERBOLIC  TANGENT 

IF  (DREAL(CTEMP).GE.  5 .D 1)  THEN 
CTANH  = (l.D0,0.D0) 

ELSE 

EXPSAV  = EXP(2.  DO*  CTEMP) 

CTANH  = (EXPSAV  - 1.D0)/(EXPSAV  + 1.D0) 

ENDIF 

DISTEL  = Rww  * CTANH  / CTEMP 
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ZrW  = rwwc+dreal(distel) 
ZiW  = dimag(distel) 


X = 1 - ZiW*Frequency*Cww 
Y = ZrW*Frequency*Cww 
DENOM  = X**2  + Y**2 
Zrww  = (ZrW*X  + ZiW*Y)/DENOM 
Ziww  = -(Ziw*X  - ZrW*  Y)/DENOM 


RETURN 

end 


c 

c 

c 

c 

c SUBROUTINE  DEBYO 1 .FOR 

c 

c This  subroutine  calculates  spectral  lineshapes  with  a Debye 
c oscillator. 

c L.  H.  Garcia-Rubio,  1991 

c Ref:  PhD  Thesis,  Me  Master  University,  1981 
c 

c w = wavenumber  (1/cm) 

c delta  = susceptibility  difference  (approx  2*epp(max)) 
c tau  = dipole  relaxation  time 

c 

c 

c 

subroutine  debye(nosc,w,tau,xip,xipp) 
real*8  zrww,ziww,rwa(  1 5),twa(  15),cwa(  1 5),freq,rwca(  1 5) 
c 

dimension  tau(4,15) 

zr=0.0 

zi=0.0 

c 


do  13  nc=l,nosc 
rwa(nc)=dble(tau(  1 , nc)) 
twa(nc)=dble(tau(2,nc)) 
rwca(nc)=dble(tau(3  ,nc)) 
cwa(nc)=dble(tau(4,nc)) 
freq=dble(w) 

call  capwarburg(freq,rwa(nc),twa(nc),cwa(nc),rwca(nc),zrww,ziww) 

zr=zr+sngl(zrww) 

zi=zi+sngl(ziww) 


13  continue 


c 

101  continue 
xip  = zr 
xipp=  zi 
c 

return 

end 

c 

c 

c 


subroutine  parfile(np) 
common/files/farr 

common/stat/ssqr,ssqi,ssq,erssq,eissq,essq,ssdatr,ssdati,ssdat 

common/osc/nosc,tau,eo,eou,eol 

common/opt/nda,nfile,nrand,ntest,nweight,mit 

common/type/fitype,comment 

common/par/bpar,upar 

common/erst/alpha,beta,gamma, delta 

common/akaike/fpe,wfpe,aic,waic 

dimension  bpar(34),upar(34),tau(4,15) 
character  farr(10)*40 
character  fitype*l,comment*75 

pi=  3.14159265 

open(22,file='parameter',status='unknown') 

write(22,*)  'Application  of  Measurement  Models  to  I.S. ' 
write(22,*) ' MEO,  PAN,  LGHR 
write(22,*) ' Dept,  of  Chemical  Engineering  ' 
write(22,*) ' University  of  Florida  ' 

write(22,*) ' Gainesville,  FI  32611  ' 

write(22,*) 1 Results  of  regression' 
write(22,*)  comment 
write(22,*) ' File(s)  =' 

do  11  i=l,nfile 

write(22,*)  farr(i) 

11  continue 

write(22,*)  'Number  of  data  sets:  ',nfile 

write(22,*)  'Number  of  measurements  per  data  set:  ',nda 

write(22,*)  'Niunber  of  Voigt  elements:  ',nosc 

write(22,*)  Residual  sum  of  squares  (real)  = ',ssqr 
write(22,*)  Residual  sum  of  squares  ( img)  = ',ssqi 
write(22,*)  Residual  sum  of  squares  (total)  = ',ssq 


write(22,*)  'Estimated  noise  level  (real)  = ',erssq 


write(22,*)  'Estimated  noise  level  ( img)  = ',eissq 
write(22,*)  'Estimated  noise  level  (total)  = ',essq 

write(22,*)  'Res  sum  of  Squares/Noise(  real)  = ',ssqr/erssq 
write(22,*)  Res  sum  of  Squares/Noise(  img)  = ',ssqi/eissq 
write(22,*)  'Res  sum  of  Squares/Noise(total)  = ',ssq/essq 


write(22,*) '%  Res  sum  of  Sq/Sq  of  Imp(  real)  = ',100.*ssqr/ssdatr 
write(22,*) '%  Res  sum  of  Sq/Sq  of  Imp(  img)  = ',100.*ssqi/ssdati 
write(22,*) '%  Res  siun  of  Sq/Sq  of  Imp(Total)  = ',100.*ssq/ssdat 

write(22,*)'fpe=',fpe,'  aic=',aic 
write(22,*)'wfpe- ,wfpe,'  waic=',waic 

if  (nweight.lt.O)  then 

write(22,"')  'fuction  values  used  for  calculatiing  the  weights' 

endif 


if  (nweight.ge.l.and.nweight.le.3)  then 

write(22,*)  'data  values  used  for  calculatiing  the  weights' 

endif 

if  (abs(nweight).eq.O)  then 

writc(22,*)  'unity  weights  used  in  regression' 

endif 


if  (abs(nweight).eq.l)  then 

write(22,*)  TJ.F.  error  structure  used  to  weight  regression' 

write(22,*) ' alpha=  '.alpha,'  beta=  '.beta,'  gamma=  '.gamma, 
1 ' delta=  '.delta 

endif 

if  (abs(nweight).eq.2)  then 

write(22,"‘)  'statistical  weighting  used  in  regression' 
endif 


if  (abs(nweight).eq.3)  then 

write(22,*)  'proportionl  weighting  used  in  regression' 
endif 


if  (abs(nweight).eq.4)  then 

write(22,',‘)  'actual  standard  deviations  used  in  regression' 
endif 


if  (abs(nweight).eq.5)  then 

write(22,*)  'weights  read  from  data,  weight' 

endif 


if  (fitype.eq."C")  then 

write(22,*) ' calculated  parameter  values  for  COMPLEX  fit ' 
endif 

if  (fitype.eq."R")  then 

write(22,*)  'calculated  parameter  values  for  REAL  fit ' 
endif 

if  (fitype.eq.'T')  then 

write(22,*)  'calculated  parameter  values  for  IMAGINARY  fit ' 
endif 


write(22,*) ' Mean  std.  deviation  95%  confidence' 
do  12  i=l,np 

dmean=(upar(i)+bpar(i))/2 
std  = (upar(i)  - bpar(i))/4. 
confd  =2.*100.*std/dmean 
write(22,1001)  dmean,std,  confd 
12  continue 

1001  format(3(el2.5,3x)) 

do  13  i=l,nosc 
write(22,*)  'Voigt  element  ',i 

write(22,*)  'R  = ’,tau(3,i),'  Ohms;  T = ',tau(4,i),'  sec' 
write(22,*) ' (C=  ',tau(4,i)/tau(3,i)/le-6,'  uF)' 


13  continue 

write(22,*)  'Solution  Resitance:  Rsol  = ',eo 

close(22) 

return 

end 


APPENDIX  B 


FORTRAN  CODE  FOR  MONTE-CARLO  SIMULATION 
This  program  uses  Monte-Carlo  simulation  to  calculate  the  95.4%  confidence 
interval  for  the  regression  performed  by  the  program  shown  Appendix  A.  This  program  is 
used  in  conjunction  with  the  program  shown  in  Appendix  A and  the  regression  algorithm 
developed  by  Prof  Garcia.  The  user  provides  the  number  of  Monte-Carlo  simulation  to  be 
performed.  The  recommended  minimum  number  of  simulations  is  5000. 
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subroutine  confidence(fitype,ndat,np) 

common/cl/refr,refi,wav 

common/par/bpar,upar 

common/modelc/zrm,zim 

common/sig/nsign,nsignsol,nsol 

common/junk/r2,gset,iset 

common/osc/nosc,tau,eo,eou,eol 

common/opt/nda,nfile,nrand,ntest,nweight,mit 

dimension  parm(34),pars(34),bpar(34),upar(34) 
dimension  datar(50000),datai(50000),eor(50000) 
dimension  parr(34, 50000), bigr(1000),bigi(1000) 
dimension  smallr(1000),smalli(1000),nsign(4,15) 
dimension  zrm(1000),zim(1000),r2(97),wav(1000) 
dimension  refr(1000),refi(1000),tau(4, 15) 
character  fitype*  1 

open(l,file='conf.real',status='unknown') 
open(2,file='conf.im',status='unknown') 
open(12, file-mean,  std', status-unknown') 
open(13,file-test.file', status-unknown') 
open(14,file='test.var',status='unknown') 

pi=3. 14159265 
nsti  = 5 
nsto  = 6 

do  111  i=l,np 


parm(i)=(upar(i)+bpar(i))/2 

pars(i)=(upar(i)-parm(i))/2 

write(12,*)  parm(i),pars(i) 

111  continue 

if  (fitype.eq.'T1)  then 

eom=(eou+eol)/2 

eos=(eom-eol)/2 

endif 


close(12) 

write(l,'(a)')  fitype 


iduml=-l 


do  21  i=l,nrand 


do  19  k=l,np 
idum=iduml 

parr(k,  i)=gasdev(idum)  *pars(k)+parm(k) 
iduml=iduml-l 
19  continue 

if  (fitype.eq.T)  then 
idum=iduml 

eor(i)=gasdev(idum)  *eos+eom 
iduml=iduml-l 

endif 

21  continue 

do  50  kk=l,nda 
c 

do  131  i=l,nrand 
w = wav(kk) 

call  debyem(nosc,w,parr,xip,xipp,i) 
c 

if  (fitype.ne."I")  then 
ep  = parr(np,i)  + xip 
endif 

if  (fitype.eq.T')  then 
ep  = eor(i)  + xip 
endif 

epp=  xipp 

datar(i)=ep 

datai(i)=epp 

if  (kk.eq.ntest)  then 
write(13,*)  i,ep,epp 
endif 

13 1 continue 


sumr=0.0 
sumi=0.0 
do  39  jj=l,nrand 
sumr=sumr+datar(jj) 
sumi=sumi+datai(ij) 
39  continue 

dmeanr=sumr/nrand 

dmeani=sumi/nrand 


sumrstd=0.0 

sumistd=0.0 


do  40  ll=l,nrand 

sumrstd=sumrstd+(datar(ll)-dmeanr)**2 
sumistd=sumistd+(datai(ll)-dmeani)**2 
40  continue 

stdr=(sumrstd/(nrand-l))**0.5 
stdi=(sumistd/(nrand- 1))*  *0. 5 

if  (kk.eq.ntest)  then 

write(14,*)  dmeanr,dmeani,stdr,stdi 

print*, dmeanr.dmeani 

print*,  stdr.stdi 

endif 


bigr(kk)=dmeanr+2*stdr 

bigi(kk)=dmeani+2*stdi 


smallr(kk)=dmeanr-2*stdr 

smalli(kk)=dmeani-2*stdi 


write(l,*)  wav(kk)  ,bigr(kk),  zrm(kk)  ,smallr(kk) 
write(2,*)  wav(kk),bigi(kk),zim(kk),smalli(kk) 


50  continue 


close(l) 

close(2) 

close(13) 

close(14) 

end 


function  gasdev(idum) 
data  iset/0/ 

common/junk/r2(97),gset,iset 
if  (iset.eq.0)  then 
1 vl=2.*ranl(idum)-l 
v2=2.  *ranl(idum)-l 


r=vl**2  + v2**2 

if(r.ge.l..or.r.eq.O.)  go  to  1 

fac=sqrt(-2.  *log(r)/r) 

gset=vl*fac 

gasdev=v2*fac 

iset=l 

else 

gasdev=gset 

iset=0 

endif 

return 

end 

function  ranl(idum) 
common/junk/r(97),gset,iset 

parameter  (ml=134456,ial=8121,icl=28411,rml=l./ml) 

parameter  (m2=243000,ia2=4561,ic2=51349,rm2=l./m2) 
parameter  (m3=259200,ia3=7141,ic3=54773) 

data  iff70/ 

if  (idum.lt.0.or.i£F.eq.0)  then 
iff=l 

ixl=mod(icl-idum,ml) 

ixl=mod(ial*ixl+icl,ml) 

ix2=mod(ixl,m2) 

ixl=mod(ial*ixl+icl,ml) 

ix3=mod(ixl,m3) 

do  11  j=l,97 

ixl=mod(ial*ixl+icl,ml) 
ix2=mod(ia2*ix2+ic2,m2) 
r(j)=(float(ixl)+float(ix2)*rm2)*rml 
11  continue 

idum=l 

endif 

ixl=mod(ial*ixl+icl,ml) 
ix2=mod(ia2‘"ix2+ic2,m2) 
ix3=mod(ia3  *ix3+ic3  ,m3 ) 
j=l+(97*ix3)/m3 
if  (j.gt.97.or.j.lt.l)  pause 

ranl=r(j) 

r(j)=(float(ixl)+float(ix2)*rm2)*rml 

return 

end 

subroutine  debyem(nosc,w,parr,xip,xipp,idat) 
common/sig/nsign,nsignsol,nsol 
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real*8  zrww,ziww,rww(15),tww(15),cww(15),rcww(15),freq 
dimension  parr(34,50000),nsign(4,15) 
c 

zr=0.0 

zi=0.0 

npar=l 

do  13  nc=l,nosc 

if  (nsign(l,nc).ne.O)  then 

rww(nc)=dble(parr(npar,idat)) 

npar=npar+l 

else 

rww(nc)=0.0 

endif 

if  (nsign(2,nc).ne.O)  then 

tww(nc)=dble(parr(npar,idat)) 

npar=npar+l 

else 

tww(nc)=0.0 

endif 

if  (nsign(3,nc).ne.O)  then 

rcww(nc)=dble(parr(npar,idat)) 

npar=npar+l 

else 

rcww(nc)=0.0 

endif 

if  (nsign(4,nc).ne.O)  then 

cww(nc)=dble(parr(npar,idat)) 

npar=npar+l 

else 

cww(nc)=0.0 

endif 

freq=dble(w) 

call  capwarburg(freq,rww(nc),tww(nc),cww(nc),rcww(nc),zrww,ziww) 

zr=zr+sngl(zrww) 

zi=zi+sngl(ziww) 

13  continue 


101 


continue 
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c 

xip  = zr 

xipp=  zi 
c 

return 

c 

end 

c 


APPENDIX  C 


FORTRAN  CODE  FOR  PROCESS  MODEL 

This  program  simulates  the  impedance  response  of  a metal  hydride  with  one 
irreversible  faradaic  reaction,  one  reversible  homogenous  reaction  and  diffusion  of 
hydrogen  into  the  bulk  of  the  metal  hydride  electrode.  A finite-difference  technique  is  used 
to  solve  the  differential  equations.  The  subroutine  ss  calculates  the  steady  state  profile  for 
the  concentration  of  hydrogen  in  the  bulk  of  the  electrode.  The  program  uses  the  BAND 
subroutine  developed  by  Newman  []  to  solve  the  equations. 
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C 

C 

C Pankaj  Agarwal  3/6/91 

C 

CC  LEVM  SUBROUTINE  KSUB 

C 

C THIS  SUBROUTINE  IS  FOR  SIMULATING  THE  IMPEDANCE  RESPONSE  OF  A METAL 
HYDRIDE 

C ONE  IRREVERSIBLE  FARAD  AIC  REACTION  ,ONE  REVERSIB  LE  HOMOGENOUS  REACTION 
AND 

C A DIFFUSION  PROCESS  HAVE  BEEN  INCORPORATED  INTO  THE  PROGRAM 
C 

C THE  SUBROUTINE  KSUB  IS  CALLED  BY  THE  MAIN  OPTIMISATION  PACKAGE  LEVM 
C 

C 'M1  IS  THE  NUMBER  OF  IMPEDANCE  DATA  PONITS 

C 'FREQ'  IS  A ONE  DIMENSIONAL  ARRAY  CONTAINING  THE  FREQUENCIES  IN  RAD/SEC 
C 'P'  IS  A ONE  DIMENSIONAL  ARRAY  CONTAINING  THE  PARAMETERS  VALUES 
C 'F  A ONE  DIMENSIONAL  ARRAY  WITH  LENGTH  2M.THE  FIRST  M POINTS  ARE  THE  REAL 
C IMPEDANCE  VALUES  AND  THE  LAST  M POINTS  ARE  THE  IMAGINARY  IMPEDANCE 
VALUES 
C 
C 
C 

C DIFF  DIFFUSION  COEFFICENT  (1/SEC)  {TO  GET  THE  REAL 

C DIFF  COEFF (CM2/SEC)  MULTIPLY  DIFF  BY  SQUARE  OF  THE  THE  ASSUMED  DIFFUSION 
C LENGTH((0.1)**2)  } 

C 

C DIFFLEN  NON  DIMENSIONAL  DIFFUSION  LENGTH  (FIXED  = 1.0) 

C 

C R1  RATE  CONST  FOR  THE  ELECTROCHEMICAL  REACTION(  1/SEC) 

C 

C R2  RATE  CONST  FOR  THE  HOMOGENOUS  REACTION  (CM3 /MOLE-SEC) 

C 

C R3  RATE  CONT  FOR  THE  REVERSIBLE  HOMOGENOUS  REACTION  (CM3 /MOLE-SEC) 

C 

C OMEGA  FREQUENCY  (RAD/SEC) 

C 

C DH  NON  DIMENSIONAL  MESH  SIZE 
C 

C TIME  TIME  AT  WHICH  THE  STEADY  STATE  PROFILE  IS  CALCULATED  (SEC) 

C 

C CAP  CAPACITANCE  (FARAD) 

C 

C CMAX  MAX  HYDROGEN  CONCENTRATION  (FIXED  = 0.09  MOLE/CM3) 

C 

C SR  SOLUTION  RESISTANCE  (OHMS) 

C 

C AREA  SURFACE  AREA  OF  THE  ELECTRODE  (CM2) 

C 

C BETA  FRACTIONAL  MAXIMUM  SURFACE  COVRERAGE  (FIXED  = 1.0) 

C 


noon  oo<oooooo  non 
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C NJ  NUMBER  OF  MESH  POINTS  ( FIXED  ) 

C 

C TAFEL  ALPHA*FARAD/R*T  ( FIXED  ) 

C 

C ALPHA  APPARENT  TRANSFER  COEFFICENT 
C 

C FARAD  FARADAY'S  CONSTANT  96487  C/EQUIV 
C 

C R UNIVERSAL  GAS  CONSTANT  8.3143  J/MOLE-DEG 
C 

C T TEMPERATURE  300K 
C 

C SCONC  SURFACE  CONCENTRATION  (DIMENSIONLESS) 

C STHETA  SURFACE  COVERAGE 
C STHETA2  1.0-THETA 
C 

C DELE  VOLTAGE  PERTURBATION  ( VOLTS  ) 

C 

C 

C 

C 

C 

C START  OF  THE  MAIN  SUBROUTINE  KSUB 
C 

SUBROUTINE  KSUB (M, FREQ, P,F) 
implicit  double  precision(a-h,o-z) 

C 

COMMON/BA/  N,NJ,A(2,2),B(2,2),C(2, 1002),D(2,5),G(2) 

1 ,X(2,2),Y(2,2),ch(2,1002) 

THIS  IS  COMMON  WITH  SUBROUTINES  bcl,bulk,bcnj,BAND  AND  MATINV 

common/par/diff,difflen,bl,rl,r2,r3, omega,  dh,time, cap, cmax 
common/temp/pp,qq,rr,beta 
common/steady/sconc,stheta,stheta2 
common/taf7tafel,dele 
dimension  p(*),f(*),freq(*) 

open(14,file='tapel4') 
open(  1 5 ,file='tape  15') 
open(16,file='tapel6') 

THESE  FILES  STORE  THE  PERTURBED  REAL  AND  IMAGINARY  CONCENTRATION 
ALUES 

AT  THREE  FREQUENCIES 

farad=96487.0 
pi=3. 141592654 
dele=0.1 
tafel=12.8 

ASSIGN  THE  PARAMETERS  VALUES  CALCULATED  FROM  THE  OPTIMIZATION  ROUTINE 
THE  LOG  OF  THE  FIRST  6 PARAMETERS  ARE  USED  FOR  OPTIMIZATION.  THE  EXTENSION 


OO  U UoO  ouuuo  uooo  ^ooouooooo 
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1 DENOTES  THAT  IT  IS  THE  LOGIO  OF  THAT  PARAMETER 

m=P(l) 

r21=p(2) 

r31=p(3) 

diffl=p(4) 

difflenl=p(5) 

capl=p(6) 

time=p(7) 

sr=p(8) 

beta=p(9) 

area=p(10) 

cmax=p(ll) 

nj=dint(p(12)) 

print*, nj 

rl=10.0**rll 

r2=10.0**r21 

r3=10.0**r31 

diff=10.**(diffl) 

difflen=10.0**(-difflenl) 

cap=10.0**(-capl) 

dh=di£flen/(nj-3) 

AS  THER  ARE  DERIVATIVE  BOUNDARY  CONDITIONS  AT  BOTH  ENDS  THERE  ARE 
TWO  IMAGE  POINTS  .HENCE  (NJ-3)  IS  USED  FOR  THE  CALCULATION  OF  THE 
MESH  SIZE 

call  ss(sconc,stheta,stheta2) 
print*,sconc,stheta,stheta2 


n=2 


START  OF  THE  IMPEDANCE  CALCULATION  LOOP  FOR  EACH  FREQUENCY 

do  111  nfreq=l,m 
omega=freq(nfreq) 
do  11  i=T,n 
do  11  k=l,nj 
c(i,k)=.00 
ch(i,k)=1.0e-15 
1 continue 

CH  COMPLEX  PERTURBED  CONCENTRATION 
C COMPLEX  CORRECTION  TO  THE  INITIAL  GUESS  OF  CH 


itmax=3 


ITMAX  MAX  NO  OF  ITERATIONS  APPLICABLE  TO  NON  LINEAR  PROBLEM 

THE  FOLLOWING  STEPS  ARE  FOR  THE  CALCULATION  OF  COMPLEX  CONCENTRATION 
USING  BAND 


no  non 
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C 

do  51  iter=l,itmax 
do  21  j=l,nj 
do  20  i=l,n 
do  20  k=l,n 
a(i,k)=0.0 
b(i,k)=0.0 
d(i,k)=0.0 

20  continue 
if(j.eq.l)  then 
call  bcl(j) 

else  if(j.eq.nj)  then 
call  bcnj(j) 
else 

call  bulk(j) 
end  if 
call  band(j) 

21  continue 


UPDATE  THE  CONCENTRATION  AT  EACH  ITERATION 

do  31  k=l,nj 
ch(  1 ,k)=ch(  1 ,k)+c(  1 ,k) 
ch(2,k)=ch(2,k)+c(2,k) 

THE  FOLLOWING  STEPS  ARE  FOR  STORING  THE  PERTURBED  CONCENTRATION 
VALUES 
C 

c if(iter.eq.itmax)  then 
c if(omega.gt.4.0e5)  then 

c zz=dh*(k-2) 
c write(14,*)  zz,ch(l,k),ch(2,k) 
c end  if 

c end  if 

c if(iter.eq.itmax)  then 

c if(omega.lt..63)  then 

c zz=dh*(k-2) 
c write(15,*)  zz,ch(l,k),ch(2,k) 
c end  if 

c end  if 

c if(iter.eq.itmax)  then 
c if(omega.lt.200.0)  then 
c if(omega.  gt.  1 50)  then 

c zz=dh*(k-2) 
c write(16,*)  zz,ch(l,k),ch(2,k) 
c end  if 

c end  if 

c end  if 

205  format(3(el8.8,3x)) 
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31  continue 
51  continue 
C 

C THE  FOLLOWING  STEPS  ARE  FOR  THE  CALCULATION  OF  THE  IMPEDANCE 
C 

C ZRT  REAL  IMPEDANCE  WITHOUT  THE  CAPACITANCE  AND  SOLUTION  RESISTANCE 
C ZIT  IMAGINARY  IMPEDANCE  WITHOUT  THE  CAPACITANCE  AND  SOLUTION 
RESISTANCE 

C ZR  REAL  IMPEDANCE  WITH  THE  CAPACITANCE  AND  SOLUTION  RESISTANCE 
C ZI  IMAGINARY  IMPEDANCE  WITH  THE  CAPACITANCE  AND  SOLUTION  RESISTANCE 
C 

thetar=(rr*ch(l,2)+qq*rr*ch(2,2)+pp*qq)/(l+qq*qq) 
thetai=(qq*rr*ch(l,2)-rr*ch(2,2)-pp)/(l+qq*qq) 
tl=dele*tafel*stheta2-thetar 
zrt=t  1 *dele/(area*rl*(t  1 **2+thetai**2))/farad 
zit=thetai  *dele/(area*  r 1 * (t  1 * * 2+thetai  * * 2))/farad 
tem= 1 -zit*  omega*  cap 
den=tem**2+(omega*cap*zrt)**2 
zr=sr+(zrt*tem+omega*cap*zrt*zit)/den 
zi=(zit*tem-zrt*zrt*omega*cap)/den 
f(nfreq)=zr 
f(nfreq+m)=zi 
111  continue 
close(14) 
close(15) 
close(16) 
return 
end 
C 
C 

C THE  SUBROUTINES  SS, SSI, NEWTON, FUND  AND  FUN  CALCULATE  THE  STEADY  STATE 
C CONCENTRATION  PROFILE  OF  HYDROGEN  IN  THE  BULK 
C 
C 

C SUBROUTINE  SS  THIS  IS  CALLED  BY  THE  MAIN  PROGRAM.  RETURNS  THE 
CALCULATED 

C VALUES  OF  SS  CONCENTRATIONS  AND  SURFACE  COVERAGE 
C 

subroutine  ss(sconc,stheta,stheta2) 
implicit  double  precision(a-h,o-z) 

common/par/di£F,difilen,bl,rl,r2,r3,  omega,  dh,time,cap,cmax 
pi=3. 141592654 

call  ssl(c) 

tl=(rl+r3*cmax*c+r2*cmax*(1.0-c)) 
stheta=(rl+cmax*r3*c)/tl 
stheta2=r2  *cmax*(  1 .0-c)/t  1 
sconc=c 

return 

end 

C 

C SUBROUTINE  SSI  MAIN  SUBROUTINE  FOR  THE  CALCULATION  OF  SS  CONCENTRATION 
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C 

C THE  BOUNDARY  CONDITION  AT  Z=0  IS  LINEARIZED  ITERATIVELY  AROUND  THE 
STEADY 

C STATE  CONCENTRATION  AND  AN  ANALYTICAL  EXPRESSION  IS  USED  FOR 
CALCULATING 

C THE  CONCENTRATION  PROFILE.  THE  PROCEDURE  IS  REPEATED  UNTIL  CONVERGENCE 
IS 

C ACHEIVED 
C 

C NMAX  MAX  NO  OF  ITERATIONS  ALLOWED  FOR  THE  ABOVE  PROCEDURE 
C 

subroutine  ssl(c) 

implicit  double  precision(a-h,o-z) 

dimension  dmu(lOOO) 

common/par/diff,difflen,bl,rl,r2,r3,omega,dh,time,cap,cmax 
pi=3. 141592654 
nmax=40 
csold=0.5 
C 

C INITIAL  GUESS  FOR  THE  SURFACE  CONCENTRATION 
C 

cs=csold 

do  113  nit=l,nmax 

terml=(rl+r3*cmax*cs+r2*cmax*(1.0-cs)) 
term2=(r3-r2)*rl*r2*cmax*(1.0-cs) 
bl=(rl*r2/terml+term2/terml**2)/(diff*0.1) 
a=(r  1 * r2/term  1 +cs*term2/term  1 * *2)/(diff*0 . 1) 

C 

C B1,A  COEFFICIENTS  OF  THE  LINEARIZED  BOUNDARY  CONDITION 
C 

c print*  ,bl 

c print*, a 

c print*, cs 

C 

C THE  FOLLOWING  STEPS  ARE  FOR  THE  CALCULATION  OF  THE  EIGENVALUES  FOR 
C THE  ANALYTICAL  SOLUTION 
C 
C 

C IF  B 1 IS  LARGE  THEN  THE  EIGENVALUES  ARE  MULTIPLES  OF  PI/2 
C 

if(bl*difllen.gt.400000.0)  then 
dmu(l)=pi/2.0 
do  11  k=2,100 
dmu(k)=dmu(k-l)+pi 
11  continue 
go  to  21 
end  if 
C 

c if(bl*difflen.lt.l.0e-15)  then 
c dmu(l)=0.0 
c do  12  k=2,100 
c dmu(k)=dmu(k-l)+pi 
c 12  continue 
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c go  to  21 
c end  if 
C 

C CALCULATE  THE  EIGENVALUES  USING  NEWTON  RAPHSON 
C 

temp=0.0 

d!omu=-pi 

himu=0.0 

guess=1.55 

if(bl*difflen.gt.  100)  guess=l. 57079 
if(bl*difflen.lt.0.001)  guess=0.05 
do  20  k=l,100 
dlomu=dlomu+pi 
himu=himu+pi 

call  newton(dlomu,himu,dmu(k),guess) 
c write(ll,*)  dmu(k) 
c print*  ,dmu(k) 
guess=dmu(k)+pi 

20  continue 

21  continue 
C 

C CALCULATE  THE  SURFACE  CONCENTRATION  USING  THE  ANALYTICAL  EXPRESSION 
C 

z=0.0 
temp=0.0 
do  30  k=l,100 

temp=temp+fun(dmu(k),time,z) 

30  continue 

c=a*(1.0-temp)/bl 

csold=cs 

cs=c 

if  (dabs((cs-csold)/csold).le.l.0e-6)  return 
113  continue 

print*, 'the  system  did  not  converge' 

return 

end 


C 

C NEWTON  RAPHSON 
C 


subroutine  newton(dlox,hix,x, guess) 
implicit  double  precision(a-h,o-z) 

pi=3. 141592654 

ncount=0 

error=1.0e-8 

numit=100 

x=guess 

do  21  i=l,numit 
call  fund(x,f,df) 
dx=t7df 


non  ooo  non 


x=x-dx 

c if((dlox-x)*(x-hix).lt.O.O)  print*,'root  out  of  bounds' 
if(dabs(dx).lt.error)  return 
21  continue 

print*,  'newton  raphson  did  not  converge' 

return 

end 


SUBROUTINE  FUND  FUNCTION  WHOSE  ROOTS  GIVE  THE  EIGENVALUES 

subroutine  fund(x,f,df) 
implicit  double  precision(a-h,o-z) 

common/par/di£f,difflen,bl,rl,r2,r3,omega,dh, time, cap, cmax 

f=x*dtan(x)-b  1 *dililen 
df=dtan(x)+x/(dcos(x)*  *2) 
return 
end 


FUNCTION  FUN  ANALYTICAL  EXPRESSION  FOR  CONCENTRATION 

function  fun(y,t,z) 

implicit  double  precision(a-h,o-z) 

common/par/diff,difflen,bl,rl,r2,r3, omega, dh,time,cap,cmax 
tl=2.0*y 

t2=-y**2*diff*t/difflen**2 
if(y.eq.0.0)  then 
en=1.0 
go  to  11 
end  if 

en=2.0*dsin(y)/(y+0.5*dsin(tl)) 

11  continue 
t3=y*(1.0-z) 

fun=en*dexp(t2)*dcos(t3) 

return 

end 


SUBROUTINE  BC1  SUBROUTINE  FOR  BAND  AT  THE  FIRST  MESH  POINT 


subroutine  bcl(j) 

implicit  double  precision(a-h,o-z) 

COMMON/BA/  N,NJ,A(2,2),B(2,2),C(2, 1002),D(2,5),G(2) 

1 ,X(2,2),Y(2,2),ch(2, 1002) 
common/temp/pp,qq,rr,beta 

common/par/di£f,difflen,bl,rl,r2,r3, omega,  dh, time, cap, cmax 


common/  steady/ scone,  stheta,  stheta2 

common/taf/tafel,dele 

pp=r  1 *dele*tafel*stheta2/(beta*omega) 

qq=(r  1 +r2  *cmax*  ( 1 . 0-sconc)+r3  *cmax  * sconc)/(beta*omega) 

rr=(r2*cmax*stheta+r3*cmax*stheta2)/(beta*omega) 

tl=(r2*(1.0-sconc)+r3*sconc)/(1.0+qq*qq) 

t2=(r2  * stheta+r3  * stheta2) 

sss=qq*rr*tl-t2 

tt=rr*tl 

uu=pp*qq*tl 

x(l,l)=-diff*0. 1/2.0 

x(l,2)=0.0 

b(l,l)=diff*0. 1/2.0 

b(l,2)=0.0 

d(l,l)=-dh*sss 

d(l,2)=dh*(tt) 

x(2,l)=0.0 

x(2,2)=-diff*0. 1/2.0 

b(2,l)=0.0 

b(2>2)=diff*0. 1/2.0 

d(2,l)=-dh*(sss) 

d(2,2)=-dh*tt 

tl=-x(l,l)*ch(l,3)-x(l,2)*ch(2,3) 

t2=-d(l,l)*ch(l,2)-d(l,2)*ch(2,2) 

t3=-b(l,l)*ch(l,l)-b(l,2)*ch(2,l)-dh*uu/qq 

g(l)=tl+t2+t3 

tl  l=-x(2,  l)*ch(l,3)-x(2,2)*ch(2,3) 

t22=-d(2,l)*ch(l,2)-d(2,2)*ch(2,2) 

t33=-b(2,l)*ch(l,l)-b(2,2)*ch(2,l)+dh*uu 

g(2)=tll+t22+t33 

return 

end 


C 

C SUBROUTINE  BULK  SUBROUTINE  FOR  BAND  IN  THE  BULK 
C 


subroutine  bulk(j) 

implicit  double  precision(a-h,o-z) 

COMMON/BA/  N,NJ,A(2,2),B(2,2),C(2, 1002),D(2,5),G(2) 

1 ,X(2,2),Y  (2,2),ch(2, 1002) 

common/par/diff,difflen,bl,rl,r2,r3,omega,dh,  time,  cap,  cmax 

a(l,l)=diff 

a(l,2)=0.0 

b(l,l)=-2*diff 

b(l,2)=-dh*dh*omega 

d(l,l)=diff 

d(l,2)=0.0 

a(2,l)=0.0 

a(2,2)=diff 


ooo  00  non 
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b(2,l)=dh*dh*omega 

b(2,2)=-2*diff 

d(2,l)=0.0 

d(2,2)=diff 

tl=-a(l,l)*ch(lj-l)-a(l,2)*ch(2j-l) 
t2=-b(l,l)*ch(l  j)-b(l,2)*ch(2j) 
t3=-d(l,l)*ch(lj+l)-d(l,2)*ch(2j+l) 
g(l)=tl+t2+t3 

tll=-a(2,l)*ch(lj-l)-a(2,2)*ch(2j-l) 

t22=-b(2,l)*ch(lj)-b(2,2)*ch(2j) 

t33=-d(2,l)*ch(lj+l)-d(2,2)*ch(2j+l) 

g(2)=tl  l+t22+t33 

return 

end 


SUBROUTINE  BCNJ  SUBROUTINE  FOR  BAND  AT  THE  LAST  MESH  POINT 


subroutine  bcnj(j) 

implicit  double  precision(a-h,o-z) 

COMMON/BA/  N,NJ,A(2,2),B(2,2),C(2, 1002),D(2,5),G(2) 

1 ,X(2,2),Y(2,2),ch(2, 1002) 

common/par/diff,di01en,bl,rl,r2,r3,omega,dh,time,cap,cmax 

do  10  i=l,2 
do  10  k=l,2 
y(i,k)=0.0 
b(i,k)=0.0 
a(i,k)=0.0 
10  continue 
a(l,l)=dh 
a(2,2)=dh 
y(l,l)=-0.5 
y(2,2)=-0.5 
b(l,l)=0.5 
b(2,2)=0.5 

tl=-y(l,l)*ch(lj-2)-y(l,2)*ch(2j-2) 
t2=-a(  1 , l)*ch(l j-  l)-a(  1 ,2)*ch(2  j-1) 
t3=-b(l,l)*ch(lj)-b(l,2)*ch(2j) 
g(l)=tl+t2+t3 

tll=-y(2,l)*ch(lj-2)-y(2,2)*ch(2j-2) 

t22=-a(2,l)*ch(lj-l)-a(2,2)*ch(2j-l) 

t33=-b(2,l)*ch(lj)-b(2,2)*ch(2j) 

g(2)=tl  l+t22+t33 

return 

end 


NEWMAN’S  BAND 
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SUBROUTINE  BAND(J) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  E(2, 3,1002) 

COMMON/BA/  N,NJ,  A(2,2),B(2,2),C(2, 1002),D(2,5),G(2) 
1 ,X(2,2),  Y (2,2), ch(2, 1002) 

101  FORMAT  (15H  DETERM=0  AT  J=,I4) 

C PRINT  J,N 
NP1  = N + 1 
IF  (J-2)  1,6,8 

1 CONTINUE 
DO  2 1=1, N 
D(I,2*N+1)  = G(I) 

DO  2 L=1,N 
LPN  = L + N 

2 D(I,LPN)  = X(I,L) 

CALL  MATINV(N,2*N+1, DETERM) 

IF  (DETERM)  4,3,4 

3 PRINT  101,  J 

4 DO  5 K=1,N 
E(K,NP1,1)  = D(K,2*N+1) 

DO  5 L=1,N 

E(K,L,1)  = - D(K,L) 

LPN  = L + N 

5 X(K,L)  = - D(K,LPN) 

RETURN 

6 DO  7 1=1, N 
DO  7 K=1,N 
DO  7 L=1,N 

7 D(I,K)  = D(I,K)  + A(I,L)*X(L,K) 

8 IF  (J-NJ)  11,9,9 

9 DO  10  1=1, N 
DO  10  L=1,N 

G(I)  = G(I)  - Y(I,L)*E(L,NPl,J-2) 

DO  10  M=1,N 

10  A(I,L)  = A(I,L)  + Y(I,M)*E(M,L,J-2) 

11  DO  12 1=1, N 
D(I,NP1)  = - G(I) 

DO  12  L=1,N 

D(I,NP1)  = D(I,NP1)  + A(I,L)*E(L,NP1,J-1) 

DO  12  K=1,N 

12  B(I,K)  = B(I,K)  + A(I,L)*E(L,K,J-1) 

CALL  MATINV(N,NP  1 , DETERM) 

IF  (DETERM)  14,13,14 

13  PRINT  101,  J 

14  DO  15  K=1,N 
DO  15  M=1,NP1 

15  E(K,M,J)  = - D(K,M) 

IF  (J-NJ)  20,16,16 

16  DO  17  K=1,N 

17  C(K,J)  = E(K,NP1,J) 

DO  18  JJ=2,NJ 

M = NJ  - JJ  + 1 


DO  18  K=1,N 
C(K,M)  = E(K,NP1,M) 

DO  18  L=1,N 

18  C(K,M)  = C(K,M)  + E(K,L,M)*C(L,M+1) 
DO  19  L=1,N 

DO  19  K=1,N 

19  C(K,1)  = C(K,1)  + X(K,L)*C(L,3) 

20  RETURN 
END 


SUBROUTINE  MATINV(N,M, DETERM) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  ID(12) 

COMMON/BA/  NTEMP,NJ,A(2,2),B(2,2),C(2,1002),D(2,5),G(2) 
1 ,X(2,2),Y(2,2)  ,ch(2,1002) 

DETERM=1.01 
DO  1 I=1,NTEMP 

I ID(I)=0 

DO  18  NN=1,NTEMP 

BMAX=1.1 

DO  6 I=1,NTEMP 

IF(ID(I).NE.O)  GO  TO  6 

BNEXT=0.0 

BTRY=0.0 

DO  5 J=1,NTEMP 

IF(ID(J).NE.O)  GO  TO  5 

IF(DABS(B(I,J)).LE.BNEXT)  GO  TO  5 

BNEXT=DABS(B(I,J)) 

IF(BNEXT.LE.BTRY)  GO  TO  5 

BNEXT=BTRY 

BTRY=DABS(B(I,J)) 

JC=J 

5 CONTINUE 

IF(BNEXT.GE.BMAX*BTRY)  GO  TO  6 

BMAX=BNEXT  /BTRY 

IROW=I 

JCOL=JC 

6 CONTINUE 
IF(ID(JC)  EQ.  0)  GO  TO  8 
DETERM=0.0 
RETURN 

8 ID(JCOL)=l 

IF(JCOL.EQ.IROW)  GO  TO  12 
DO  10  J=1,NTEMP 
SAVE=B(IROW,J) 

B(IROW,J)=B(JCOL,J) 

10  B(JCOL,J)=SAVE 
DO  11  K=1,M 
SAVE=D(IROW,K) 

D(IROW,K)=D(JCOL,K) 

II  D(JCOL,K)=SAVE 
12  F=  1.0/B(JCOL, JCOL) 

DO  13  J=1,NTEMP 
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13  B(JCOL,J)=B(JCOL,J)*F 
DO  14  K=1,M 

14  D(JCOL,K)=D(JCOL,K)  *F 
DO  18  I=1,NTEMP 
IF(I.EQ.JCOL)  GO  TO  18 
F=B(I,JCOL) 

DO  16  J=1,NTEMP 

16  B(I,J)=B(I,J)-F*B(JCOL,J) 
DO  17  K=1,M 

17  D(I,K)=D(I,K)-F*D(JCOL,K) 

18  CONTINUE 
RETURN 
END 
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